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Due to the development of efficient algorithms and the im- 
provement of computer power it is now possible to map out 
potential energy surfaces (PES) of reactions at surfaces in 
great detail. This achievement has been accompanied by an 
increased effort in the dynamical simulation of processes on 
surfaces. The paradigm for simple reactions at surfaces - 
the dissociation of hydrogen on metal surfaces - can now be 
treated fully quantum dynamically in the molecular degrees 
of freedom from first principles, i.e., without invoking any 
adjustable parameters. This relatively new field of ab ini- 
tio dynamics simulations of reactions at surfaces will be re- 
viewed. Mainly the dissociation of hydrogen on clean and ad- 
sorbate covered metal surfaces and on semiconductor surfaces 
will be discussed. In addition, the ab initio molecular dynam- 
ics treatment of reactions of hydrogen atoms with hydrogen- 
passivated semiconductor surfaces and recent achievements 
in the ab initio description of laser-induced desorption and 
further developments will be addressed. 
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Understanding reactions on surfaces plays an impor- 
tant role in a wide range of technologically relevant ap- 
plications. Among those are the heterogenous catalysis 
- the majority of reactions in the chemical industry em- 
ploy catalysts; crystal growth, which determines, e.g., the 
quality of semiconductor devices; corrosion and lubrica- 
tion, which influences the durability of mechanical sys- 
tems; or hydrogen storage in metals, just to mention a 
few. The reactions involved in these processes are of- 
ten too complicated to be studied in detail as a whole. 
Therefore in surface science one tries to understand re- 
action mechanisms by breaking them up into simpler 
steps which are then studied under well-defined condi- 
tions |,|. 

Studies of chemical reaction dynamics are well estab- 
lished in the gas phase. Reactions on surfaces differ from 
gas-phase reactions in two fundamental aspects: first, 
the presence of the surface changes the symmetry signif- 
icantly, and second, the substrate on which the reaction 
occurs represents in principle a system with an infinite 
number of degrees of freedom which acts as a heat bath 
that leads to dissipation and thermal fluctuations. As for 
the symmetry, while for example for the description of the 
dissociation of a diatomic molecule in the gas phase only 
the interatomic distance has to be taken into account 
(all other degrees of freedom can be separated), in front 
of a surface all six degrees of freedom of the diatomic 
molecule have to be considered explicitly since the sur- 
face breaks the symmetry with respect to the gas phase. 
On the other hand, a crystalline surface introduces new 
symmetries, for example the periodicity of the surface. 

Dissociative adsorption processes are of particular im- 
portance for reactions on surfaces since they constitute 
the first step in heterogenous catalysis; furthermore they 
are often the rate-limiting process, for example in ammo- 
nia synthesis or CO oxidation. The term "heterogenous" 
refers to the fact that the reactants (usually a gas) and 
the catalyst (usually a solid substrate) are not in the 
same phase state. For the atomic or molecular adsorp- 
tion process, energy dissipation to the substrate is nec- 
essary for the particles to be trapped into an adsorption 
well on the surface, otherwise they would be scattered 
back into the gas phase. In the case of dissociative ad- 
sorption, however, there is another channel for energy 
transfer, which is the conversion of the kinetic and inter- 
nal energy of the molecule into translational energy of the 
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atomic fragments on the surface relative to each other. In 
addition, for light molecules like hydrogen there is usu- 
ally only a small energy transfer to substrate phonons 
due to the large mass mismatch between the molecule 
and the substrate atoms. If furthermore no substantial 
surface rearrangement upon adsorption occurs - which is 
often fulfilled for densely packed metal surfaces - then the 
dissociative adsorption dynamics can be described using 
low-dimensional potential energy surfaces (PES). 

Still the interaction of a diatomic molecule with a well- 
defined fixed substrate involves the six degrees of freedom 
of the molecule. Up to recently the interaction dynam- 
ics in particular of hydrogen had been described in low- 
dimensional quantum studies |p|-^7|; higher-dimensional 
studies could only be performed by classical trajectory 
calculations [^8|-^0) or mixed quantum-classical methods 
fjl| [m| on model potentials. This was due to the lack 
of reliable potential energy surfaces and the large com- 
putational effort for high-dimensional quantum dynam- 
ics studies (fo r reviews on these dynamical studies see, 
e.g., Refs. P5|-|3^|). This situation has changed signif- 
icantly within the last five years caused by the devel- 
opment of efficient algorithms and the increase in com- 
puter power. For the paradigm of simple reactions on 
surfaces - the dissociation of hydrogen on metal surfaces 
- detailed potential energy surfaces obtained by density 
functional theory calculations are now available |39] p9| . 
This has increased the motivation to perform dynami- 
cal studies in which all hydrogen degrees of freedom are 
treated explicitly, and indeed the first six-dimensional 
quantum dynamical studies of hydrogen dissociation on 
metal surfaces have now been performed [j5^ . 

In this review I will give a overview of this rather new 
field of ab initio dynamics calculations of reactions on 
surfaces. Usually these calculations require three inde- 
pendent steps: 

1. Determination of the ab initio PES by first- 
principles total-energy calculations, 

2. a fit of the total energies to an analytical or numeri- 
cal continuous representation which serves as an in- 
terpolation between the actually calculated points, 

3. a dynamical calculation on this representation of 
the ab initio PES that includes all relevant degrees 
of freedom. 

Since this type of calculations is indeed derived from first 
principles with no adjustable parameters, I will refer to 
them as ll ab initio dynamics calculations". However, this 
term should not conceal the still approximative nature of 
these dynamics simulations. In all three steps approxima- 
tions are involved which will be discussed in this review. 
For example, exchange and correlation effects of the elec- 
trons can not be treated exactly in total-energy calcula- 
tions using density functional theory. The description of 
the surface by either a supercell or a cluster method is an 
approximation. Furthermore, the fitting of the ab initio 



energies to a continuous representation is a highly non- 
trivial task. If just classical trajectories are determined, 
quantum effects in the motion of the nuclei are neglected 
which can be important, in particular for the dynamics of 
hydrogen. And also the question, what are the relevant 
degrees of freedom in the dissociation process, cannot be 
answered unambiguously and depends, for example, on 
the properties that one is interested in. Hence in this 
review also investigations, in which not all molecular de- 
grees of freedom are treated explicitly, but which are still 
derived from first principles, are included. 

For classical ab initio dynamics simulations, which will 
in the following be referred to by the term 11 ab initio 
molecular dynamics" , the three steps mentioned above 
can actually be combined since most first-principles total- 
energy schemes also determine the gradients of the po- 
tential via the Hellmann-Feynman theorem Q,^5| . Thus 
the classical equations of motion can be directly solved. 
For each step of the numerical integration of the equa- 
tions of motion the forces are determined by a new total- 
energy calculation. Since this is rather time-consuming, 
only a small number (< 100) of trajectories can be calcu- 
lated in this way. Therefore this method does not allow 
the determination of reaction probabilities which usually 
require at least thousands of trajectories due to the sta- 
tistical nature of the reactive events |5l| . 

On the other hand, ab initio molecular dynamics sim- 
ulations with the determination of the forces "on the 
fly" do not require any fitting. This makes them very 
flexible. On semiconductor surfaces, the substrate rear- 
rangement upon adsorption cannot be neglected; indeed 
it plays a crucial role for the adsorption and desorption 
mechanism |^J56|] . An analytical fit of an ab initio PES 
including several substrate degrees of freedom has not 
been performed yet because of its complexity. Therefore, 
for the investigation of the dissociative adsorption and 
associative desorption process on semiconductor surfaces 
the "traditional" ab initio molecular dynamics method 
with the determination of the forces "on the fly" has still 
been used [p7j-p0| . Its application is reasonable if the in- 
formation obtained from a small number of trajectories 
is sufficient to gain insight into a particular process. 

High-dimensional ab initio dynamics calculations will 
not make low-dimensional simulations on model poten- 
tials obsolete. These model calculations have provided 
us with the framework to interprete the dynamics. High- 
dimensional dynamics is often too complicated to be fol- 
lowed in detail. But this review will show that high- 
dimensional ab initio dynamics calculations have ad- 
vanced our understanding of the dissociation dynamics 
tremendously and sometimes even caused the modifica- 
tion of established concepts. First of all, the ab initio 
calculation of potential energy surfaces has confronted 
the gas-surface dynamics community with some features 
of these PESs which had not been expected and which 
challenged new interpretations and calculations. And 
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secondly, some phenomema in the reaction dynamics on 
surfaces only occur in simulations if a sufficiently large 
number of degrees of freedom is included. Thus high- 
dimensional simulations not only lead to progress in the 
quantitative, but also in the qualitative understanding of 
processes on surfaces. In particular for the hydrogen dis- 
sociation on transition metals surfaces ab initio dynam- 
ics calculations have proven their power. They demon- 
strated, e.g., that the initial decrease of the sticking prob- 
ability with kinetic energy often found experimentally in 
these systems is not due to a precursor mechanism, as 
was the common believe, but is caused by an hitherto un- 
derestimated dynamical mechanism, namely the steering 
effect HIOJ. 

This review is devoted to a large part to the dissocia- 
tion of hydrogen on surfaces. From the theoretical point 
of view, hydrogen is the simplest molecule. This does not 
neccessarily keep the computational effort small, but it 
makes the theoretical description controllable. From the 
experimental point of view, hydrogen molecules can be 
used without severe problems in molecular beam appara- 
tuses, and they can be detected state-specifically by laser 
techniques. Therefore hydrogen is the ideal candidate for 
a close collaboration between theory and experiment in 
order to explore the dynamics of the dissociation process 
on surfaces. The general concepts found in these studies 
will also be applicable to heavier and more complicated 
molecules. Due to its light mass, the study of the hydro- 
gen dissociation dynamics also allows to address some 
fundamental concepts in physics like the importance of 
quantum effects. Still the hydrogen interaction with sur- 
face is technologically relevant, for example in the pas- 
sivation or growth of semiconductor devices. All these 
facts make the study of the hydrogen dissociation dy- 
namics on surfaces an exciting research subject. 

The electronically adiabatic ab initio dynamics of ad- 
sorption and desorption has recently been reviewed from 
a more quantum chemical point of view [|60| . In the 
meantime, an ab initio molecular dynamics study of the 
reaction of hydrogen atoms with a hydrogen-passivated 
semiconductor surface has been performed [B2| in which 
pick up reactions of the Eley-Rideal type were studied. 
Furthermore, detailed ab initio PES studies of reactions 
on surfaces involving heavier molecules like O2 |^3| or 
CO2 nave been performed, and probably soon dy- 

namics calculations using these ab initio potentials will 
be carried out. Even the laser-induced desorption which 
involves electronically excited states has currently been 
addressed from first principles, namely the laser-induced 
desorption of NO form NiO(lOO) US which wiU also 
be briefly reviewed at the end of this article. This short 
collection shows that ab initio dynamics studies of reac- 
tions on surface is a very active and growing field. 

After this introduction, in the next section the theo- 
retical concepts necessary in order to perform ab initio 
dynamics calculations are briefly introduced. Then the 



hydrogen dissociation on metal surfaces is addressed, first 
on simple and noble metals where the dissociation is usu- 
ally hindered by a substantial barrier, and then on tran- 
sition metals which often have non-activated pathways 
to dissociative adsorption. In this context also the poi- 
soning of the dissociation by an adsorbate is discussed. 
The interaction of hydrogen with semiconductor surfaces 
is treated in the following section. Then the first ab ini- 
tio dynamics calculation involving electronically excited 
states is discussed. The review ends with an outlook and 
some concluding remarks. 

II. THEORETICAL BACKGROUND 

In this section the theoretical methods to describe the 
molecule-surface interaction are briefly described. It is 
not intended to give a complete overview. I will rather in- 
troduce the main concepts and shortly comment on them. 

A. Born-Oppenheimer Approximation 

The Schrodinger equation describing the interaction of 
molecules with surfaces has the general form 

H |*({R m ,r„})) = £ ^({B^.rn})). (1) 

In this equation R m are the ionic coordinates and r„ 
the electronic coordinates. It is well-known that a com- 
plete analytical solution of the Schrodinger equation tak- 
ing into account both ionic and electronic degrees of free- 
dom is not possible except for simple cases. One common 
approach is to assume that - due to the large mass differ- 
ence between electrons and the nuclei - the electrons fol- 
low the motion of the nuclei adiabatically. This is the fa- 
mous Born-Oppenheimer approximation frofl . And since 
it is only an approximation, its validity has to be checked 
carefully. 

In gas-surface scattering electronically non-adiabatic 
processes are indeed occuring. They are directly ob- 
servable as chemiluminescence and exo-electron emis- 
sion (see, e.g., the recent review of T. Greber fri|). 
However, the proper treatment of electronically non- 
adiababic effects in gas-surface dynamics is rather com- 
plex. One problem that goes with the complexity of the 
non-adiabaticity is that it is not easy to judge whether 
the dynamics in any particular system is indeed electron- 
ically adiabatic or not. The usual argument for the valid- 
ity of the Born-Oppenheimer approximation is the small- 
ness of the atomic velocities as compared to the electronic 
velocities. If in addition the adiabatic energy levels are 
well-separated, electronic transitions are negligible (see, 
e.g., Ref. For metals the situation is more com- 

plicated due to the quasi-continuum of electronic states. 
Since the effective potentials and the coupling between 
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the electronic states can still not be computed rigorously, 
one is left with more or less hand- waving arguments. 
Among them are the following: (i) At metal surfaces 
electronic excitations are very short-lived. Hence elec- 
tronic excitations are effectively quenched, (ii) Molecular 
electronic levels become rather broad upon the interac- 
tion with surfaces. Broad levels correspond to short life- 
times of excited states and again, this leads to an effective 
quenching of electronic excitations. 

On semiconductor and insulator surfaces the situation 
is different due to the band gap which makes the treat- 
ment of electronically excited states more tractable. At 
the end of this review 1 will present the first descrip- 
tion of electronically non-adiabatic processes in molecu- 
lar desorption from oxide surfaces on an ab initio basis. 
For metal surfaces, where no band gap exists, a more 
practical approach has still to be applied which is just 
to perform electronically adiabatic dynamical studies of 
reactions on surfaces as well and detailed as possible and 
to compare the results of these high-quality calculations 
with experiment. As long as the consequences of these 
simulations are in agreement with experiment, there is 
apparently no need for the involvement of electronic ex- 
citations in the described processes. For the time being, 
we just follow this practical approach and assume that 
the Born-Oppenheimer approximation is justified. We 
can then write down the electronic Hamiltonian in which 
the coordinates of the nuclei just enter as parameters. 
This Hamiltonian has the form 

H el ({R m }) |V({r„») = E ({R m }) |V({r„})>. (2) 

The many-electron ground state energy fi ({R m }) then 
defines the potential for the motion of the nuclei. Once 
it is obtained, it can be plugged into the Schrodingcr 
equation for the nuclei, 

(E^ v r. + £o({R m })) $({R™})= ^W, 

(3) 

where £ is now the energy relevant for the dynamics of 
the nuclei, or, alternatively, it can be used to solve the 
classical equations of motion, 
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For extended systems the most efficient approach 
to determine the many-electron ground state energy 
£o({Rm}) from first principles is density functional the- 
ory [ |73| , [r4| in combination with the supercell concept. In 
this approach the surface is modelled by a periodic ar- 
ray of slabs separated by vacuum regions, and the wave 
functions are expanded in a suitable basis which is often 
based on plane waves or augmented plane waves. Almost 



all of the ab initio potentials discussed later in this re- 
view are determined in this way, but real-space [j75 76 
and Green function methods J77|-[7S|] have been proposed 
as well. And also quantum chemical methods are used in 
which the infinite substrate is modelled by a finite cluster 
(see, e.g., the recent review of Whitten and Yang JSOfl). 
Ab initio total energy calculations are briefly discussed 
in the next section. 



B. Ab initio total energy calculations 

In this review I will focus on the reaction dynamics on 
ab initio potential energy surfaces, not on the determina- 
tion of the PESs. Still I will introduce some basics about 
Density Functional Theory (DFT). For further details I 
refer to more extensive treatments like, e.g., Refs. fjf| , jS2f . 
DFT is based upon the Hohenberg-Kohn theorem |73[| 
which states that the ground-state total energy of a sys- 
tem of interacting electrons -Etot can be obtained by min- 
imizing an energy functional E[n], 



min E\n] 



min(T[n] + U[n] + E xc [n}). (5) 



T[n] and U[n] are the functionals of the non- interacting 
many electron kinetic and electrostatic energy, respec- 
tively. All quantum mechanical many-body effects are 
contained in the so far unknown exchange-correlation 
functional E xc [n]. The electron density n(r) which min- 
imizes the total energy can be found by solving self- 
consistently the Kohn-Sham equations M4i 



-ti 2 
2m 



V 2 + V cs (r) + V xc (r) 



V>i(r) = £i"0i(r), 
(6) 



where T is the single-particle kinetic energy operator 
and V es (r) the electrostatic potential. The exchange- 
correlation potential V xc (r) is the functional derivative 
of the exchange-correlation functional E xc [n] 



F xc (r) 



SE X 
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(7) 



The exchange-correlation functional E xc [n] can be writ- 
ten as 



E x 



d 3 r n(r) e xc [n](r), 



(8) 



where e xc [n] is the exchange-correlation energy per par- 
ticle. This exchange-correlation energy is not known ex- 
actly except for some special cases like constant electron 
density. In a wide range of bulk and surface problems the 
so-called local density approximation (LDA), in which 
e xc [n] is replaced by the exchange-correlation energy for 
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the homogeneous electron gas, has been surprisingly suc- 
cessful |Q. For dissociation barriers on surfaces, how- 
ever, LDA is seriously in error. In the generalized gradi- 
ent approximation (GGA) also the gradient of the density 
is included in the exchange-correlation functional. Dif- 
ferent forms for the GGA functional have been proposed 
(see, e.g., Refs. |8^-|8(|). The use of GGA functional 
leads to a significant improvement in the accuracy of cal- 
culated barrier heights p9||85| |. Still the validity of the 
GGA is strongly debated, in particular in the H2/Si(100) 
system p7| . The problem in the development of more ac- 
curate exchange-correlation functionals is that they still 
represent in principle an uncontrolled approximation, i.e., 
there is no systematic way of improving the functionals 
since there is no expansion in some controllable param- 
eter (it should be noted here that even if an expansion 
in some controllable parameter exists, this expansion can 
still be problematic). Basically the success justifies the 
choice of some particular functional. Hence dynamical 
simulations on ab initio potentials also serve the purpose 
to check the accuracy of the chosen exchange-correlation 
functional. 

The supercell approach allows to transfer DFT algo- 
rithms, that have been very successfully applied to the 
determination of bulk properties, to the description of 
surface problems. Especially the seminal paper by Car 
and Parrinello has drawn a lot of attention. Due 
to the variational principle of density functional theory 
the determination of the electronic ground state can be 
regarded as a global optimization scheme thus avoiding 
the explicit diagonalization of huge matrices. Such an ap- 
proach was first proposed by Bendt and Zunger ]89| ] . Car 
and Parrinello demonstrated how powerful this method 
can be by combining DFT and molecular dynamics. 

In the quantum chemical approach the infinite sub- 
strate is represented by a finite cluster. DFT methods 
have also been used to determine the total energies of the 
cluster , but traditionally quantum chemical meth- 

ods based on the wave function dominate, i.e. Hartree- 
Fock methods without or with the inclusion of the con- 
figuration interaction (CI). While supercell calculations 
often describe relatively high-coverage situations due to 
the repeated surface unit cell, cluster methods represent 
low-coverage situations. However, they suffer from the 
fact that the infinite substrate is often not appropriately 
modelled by a small cluster and that the total energies 
are often not converged as a function of the cluster size 
fl30[ . For the treatment of electronic excitations, DFT 
methods could in principle also be used because due to 
the Hohenberg-Kohn theorem the electronic ground state 
density uniquely defines the external potential and 
thereby also all electronically excited states. However, in 
practice DFT methods usually only determine the total 
energy for the electronic ground state of a certain config- 
uration (there are exceptions, see for example Ref. [p2|). 
To describe electronically excited states wave-function 



based methods are still unavoidable which necessitates 
the use of a finite cluster that can, however, be embed- 
ded for example in a field of point charges pq,p9] . 



C. Parametrization of the ab initio potential 

One serious problem arises for the use of ab initio po- 
tential energies in particular in quantum dynamics sim- 
ulations. To solve the Schrodingcr equation, one needs 
in general a continuous description of the potential since 
the wave functions are delocalized. The ab initio calcu- 
lations, however, just provide total energies for discrete 
configurations of the nuclei. In order to obtain a contin- 
uous description, the ab initio energies have to be fitted 
to an analytical or numerical continuous representation 
of the potential energy surface. This is a highly non- 
trivial task. On the one hand the representation should 
be flexible enough to accurately reproduce the ab initio 
input data, on the other hand it should have a limited 
number of parameters so that it is still controllable. Ide- 
ally a good parametrisation should not only accurately 
interpolate between the actually calculated points, but 
it should also give a reliable extrapolation to regions of 
the potential energy surface that have actually not been 
determined by the ab initio calculations. 

The explicit form of the chosen analytical or numeri- 
cal representation of the ab initio potential varies from 
application to application. Often the choice is dictated 
by the dynamics algorithm in which the representation is 
used. The applications have almost entirely been devoted 
to the interaction of a diatomic molecule with the sur- 
face. The angular orientation of the molecule has usually 
been expanded in spherical harmonics and the center-of- 
mass coordinates parallel to the surface in a Fourier series 

lH,|o]]6l]Il. For the PES in the plane of the molecu - 
lar distance from the surface and the interatomic sepa- 
ration a representation in reaction path coordinate has 
been employed 93 94j, but also two-body poten- 

tials have been used |}47| . Before detailed ab initio po- 
tentials became available, the LEPS form was often used 
to construct a global PES fl35| . This parametrization 
contains only a small number of adjustable parameters 
which made it so attractive for model calculations, but 
which makes it at the same time relatively unflexible. A 
modified LEPS potential has still been successfully used 
to fit an ab initio PES of the interaction of atomic hy- 
drogen with the the hydrogenated Si(100) surface (62j. 

As an alternative approach to fit ab initio energies, a 
genetic programming scheme has recently been proposed, 
in which both the best functional form and the best set of 
parameters are searched for |95) . This method has so far 
only been used for three-dimensional potentials so that 
its capability still has to be proven for higher-dimensional 
problems. 
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Ab initio total energies are often mainly determined 
at high-symmetry points of the surface in order to re- 
duce the computational cost. It is true that these high- 
symmetry points usually reflect the extrema in the PES. 
However, due to this limitation the fitted continuous PES 
can only contain terms that correspond to these high- 
symmetry situations. On the one hand this often saves 
computer time also in the quantum dynamics because 
certain additional selection rules are introduced which re- 
duces the necessary basis set 47,53]. On the other hand, 
of course this represents an approximation. The question, 
how serious the neglect of terms with lower symmetry is, 
remains open until these terms have been determined and 
included in actual dynamical calculations. 

If more than just the molecular degrees of freedom 
should be considered in a parametrization of an ab initio 
PES, analytical forms become very complicated and cum- 
bersome. As an alternative, the interpolation of ab initio 
points by a neural network has been proposed [p6|-p8|. 
Neural networks can fit, in principle, any real- valued, 
continuous function to any desired accuracy. They re- 
quire no assumptions of the functional form of the un- 
derlying problem. On the other hand, there is no physi- 
cal insight that is used as an input in this parametriza- 
tion. Hence the parameters of the neural network do not 
reflect any physical or chemical property. Another ap- 
proach is to fit the parameters of a tight-binding formal- 
ism to ab initio energies p^ |. A tight-binding method is 
more time-consuming than an analytical representation 
or a neural network since it requires the diagonalization 
of a matrix. However, due to the fact that the quan- 
tum mechanical nature of bonding is taken into account 



[100 1 tight-binding schemes need a smaller number of ab 



initio input p oints to perform a good interpolation and 
extrapolation [101]. In add ition, their parameters, the 
Slater-Koster integrals [ 102 1 , have a physical meaning. It 
remains to be seen how useful the different fitting meth- 
ods will be for performing ab initio dynamics simulations. 

It has been argued that any fitting scheme using ana- 
lytic potential functions introduces a possible bias into 
the parametrization fjof . In order to avoid this one 
should use ab initio molecular dynamics methods that 
calculate the forces "on the fly" because this requires 
no fitting. However, approaching a problem apparently 
unbiased also means not to take advantage of any previ- 
ous knowledge. On the contrary, using as much experi- 
ence and knowledge as possible should usually make any 
method more efficient fl95[ . 

An important issue is to judge the quality of the fit to 
an ab initio PES. Usually the root mean squared (RMS) 
error between fit and input data is used as a measure 
of the quality of a fit. If this error is zero, then every- 
thing is fine. However, normally this error is larger than 
zero. The systematic error of the ab initio energies is 
usually estimated to be of the order of 0.1 eV. Often it is 
said that the RMS error of the fit should be of the same 



order. But the dynamics of molecular dissociation at sur- 
faces can be dramatically different depending on whether 
there is a barrier for dissociation of height 0.1 eV or not 



[ 103 1 . Hence for certain regions of the PES the error has 
to be much less than 0.1 eV, while for other regions even 
an error of 0.5 eV might not influence the dynamics sig- 
nificantly. Another example occurs in a reaction path 
parametrization. If the curvature in the parametrization 
is off by a few percent, the energetic distribution of bar- 
rier heights is not changed and the dynamical properties 
are usually not altered significantly. However, the lo- 
cation of the barriers is changed and consequently the 
RMS error can become rather large. Hence one has to be 
cautious by just using the RMS error as a quality check 
of the fit. Unfortunately there is no other simple error 
function for the assessment of the quality of a fit. If it is 
possible, one should perform a dynamical check. Obvi- 
ously, if the dynamical properties calculated on a fitted 
PES agree with the ones calculated on the original PES, 
the quality of the fit should be sufficient. 



D. Quantum dynamics 

There are two ways to determine quantum mechani- 
cal reaction probabilities: by solving the time-dependent 
or the time-independent Schr odinger equation. Both ap- 
proaches are equivalent [ 104 1 and should give the same 
results. The chosen method depends on its applicability, 
but apparently it is often also a matter of training and 
personal taste. 

In the most common time-independent formulation, 
the concept of defining one specific reaction path coordi- 
nate is crucial. Starting from the Schrodinger equation 



(H - E) * = 0, 



(9) 



one chooses one specific reaction path coordinate s and 
separates the kinetic energy operator in this coordinate 



2(1. 



H -E)V 



0. 



(10) 



Here H is the original Hamiltonian except for the kinetic 
energy operator in the reaction path coordinate. Usually 
the use of curvelinear reaction path coordinates results 
in a more complicated expression for the kinetic energy 
operator involving cross terms [105-107], but for the sake 

As the next 
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of clarity I have neglected this in Eq. 
step one expands the wave function in the coordinates 
perpendicular to the reaction path coordinate in some 
suitable set of basis functions, 



* =*(«,...) = 



(ii) 



Here n is a multi-index, and the expansion coefficients 
tpn(s) are assumed to be a function of the reaction path 



G 



coordinate. Now we insert the expansion of VP in Eq. |i0| 
and multiply the Schrodinger equation by (m|, which 
corresponds to performing a multi-dimensional integral. 
Since the basis functions |n) are assumed to be indepen- 
dent of s, we end up with the so-called coupled-channel 
equations, 



^ .(— d* - E) 5 m ,n+ {m\H\n) \ ^„(s) = 0. 

(12) 

Instead of a partial differential equation - the original 
time-independent Schrodinger equation Eq. ^| - we now 
have a set of coupled ordinary differential equation. Still 
a straightforward numerical integration of the coupled- 
channel equations leads to instabilities, except for in sim- 
ple cases, due to exponentially increasing so-called closed 
channels. Recently a very stable and ef ficient co upled- 
channel algorithm has been introduced | 10§j - 111]. The 
main idea underlying this particular algorithm will be 
briefly sketched in the following. 

For the solution VP defined in Eq. [H], which represents 
a vector in the space of the basis functions, the initial 
conditions are not specified. This function can also be 
considered as a matrix 



■III J 



(13) 



where the index I labels a solution of the Schrodinger 
equation with an incident plane wave of amplitude one 
in channel I and zero in all other channels. Formally one 
can then write the solution of the Schrodinger equation 
for a scattering problem in a matrix notation as 



*(s -> +oo) = e~ lqs - 
*(s -> -oo) = e- iqs t. 



(14) 



Here q — qmfim.n is a diagonal matrix, r and t are the 
reflection and transmission matrix, respectively. Now one 
makes the following ansatz for the wave function, 



9(s) = (l-p(s)) 



1 



r(s) 



t. 



(15) 



Equation [l5] defines the local reflection matrix p(s) 
(LORE) and the inverse local transmission matrix r(s) 
(INTRA). The boundary values for these matrices are 
(except for phase factors which, however, do not affect 
the transition probabilities): 



(p(s);r(s)) = 



(r-t) 

(0;i) 



+oo 
— oo 



(16) 



From the Schrodinger equation first order differential 
equations for both matrices can be derived |108| 
which can be solved by starting from the known initial 
values at s — > — oo; at s — > +oo one then obtains the 



physical reflection and transmission matrices. Thus the 
numerically unstable boundary value problem has been 
transformed to a stable initial value problem. 

In the time-dependent or wave-packet formulation, the 
solution of the time-dependent Schrodinger equation 



i% *(R,t) = H #(R,t) 
can formally be written as 

#(R,t) = e ~ iHt/n *(R,t = 0), 



(17) 



(18) 



if the potential is time-independent. The most com- 
mon methods to represent the time-evolution opera- 
tor exip(—iHt/h) in the g as-surfac e dynamics commu- 
nity are the split-operator 1 112] , 113(1 an d the Chebychev 



114 1 methods. In the split-operator method, the time- 



evolution operator for small time steps At is written as 



-iHAt/h 



-iKAt/2h e -iVAt/h e ~iKAt/2h _|_ Q(/±f3^ 

(19) 



where K is the kinetic energy operator and V the poten- 
tial term. In the Chebyshev method, the time-evolution 
operator is expanded as 



-iHAt/h 



a, (At) Tj(H), 



(20) 



where the Tj are Chebyshev polynomials and H is 
the Hamiltonian rescaled to have eigenvalues in the 
range (-1,1). Both propagation schemes use the fact that 
the kinetic energy operator is diagonal in k-space and 
the potential in real-space. The wave function and the 
potential are represented on a numerical grid, and the 
switching between the k-space and real-space representa- 
tions is effic iently done by Fast Fourier Transformations 
(FFT) fn|]. 

In the last years the time-dependent wave-packet 
methods have been much more fashionable than time- 
independent schemes in the gas-surface dynamics com- 
munity. Many different groups have used wave-packet 
codes to study the dissociation dynamics on surfaces, 
up to recently almost entirely on low-dimensional model 
potentials j| |,p^p^|Jll|Jl^ p6l^ll . It has been argued 
that wave-packet methods avoid "the problem of exces- 
sively many channels" |3j|. But actually these methods 
suffer from the fact that the wave function and the po- 
tential are represented on a grid which leads to memory 
problems in the implementation. In each dimension be- 
tween 16 and several hundred grid-points are used. In 
time-independent methods, on the contrary, it is for ex- 
ample sufficient to expand the hydrogen wave function 
in the interatomic distance in three or four eigenfunc- 
tions of a harmonic oscillator due to the large vibra- 
tional energy quantum. One might say that the hy- 
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drogen vibration in time-independent methods is repre- 
sented by three or four points instead of more than six- 
teen in time-dependent methods. In addition, while in 
the LORE - INTRA scheme p8| - |lTl| the matrices are 



successively determined along the reaction path and do 
not have to be stored, in wave-packet methods the wave 
function has always to be stored everywhere on the grid. 
The saving in storage requirements is one of the reasons 
that the first six-dimensional quantum dynamical treat- 
ment of hydrogen dissociation on surfaces was indeed a 
time- independent one using the LORE— INTRA scheme 



E. Classical dynamics 

Once an analytical continuous representation of an ab 
initio PES is available, it is also possible to determine 
the gradients of the potential analytically. This allows 
one to perform ab initio molecular dynamics studies by 
integrating the classical equations of motion, 



d 2 



d 



(21) 



The solution of the equations of motion can be obtained 
by standard numerical integration schemes like Runge- 
Kutta, Burlisch-Stoer or predictor-corrector methods 
(see, e.g., Ref. |115| ]). Sticking corresponds to a pro- 
cess in which statistically distributed particles hit the 
surface. This means that the determination of classical 
sticking probabilities requires an average over typically 
thousands of trajectories. The initial conditions can be 
chosen either by some Monte-Carlo sampling or by an 
equidistant sampling. 

Most ab initio total-energy programmes based on 
plane-wave expansions also evaluate the gradients of the 
potential via the Hellmann-Feynman theorem p4| , |55| |. 
With these forces the classical equations of motion can 
directly be solved. This "traditional" ab initio molecular 
dynamics scheme with the determination of the forces 
on the fly has the advantage that it does not require 
the fitting of the ab initio PES to any analytical or nu- 
merical representation. This is especially advantageous 
if surface degrees of freedom play an important role in 
the dissociation process. On the other hand, every step 
of the numerical integration of the equations of motion 
requires a new ab initio total-energy calculation which is 
still rather time-consuming. Hence the number of trajec- 
tories obtainable in such a "during the journey" ab ini- 
tio molecular dynamics simulation is limited to numbers 
well below 100 f37|-|60|]. Such a simulation only makes 
sense if one can extract meaningful information from a 
low number of events. For the simulation of adsorption 
events this is usually not the case, but for desorption, 
where the particles originate from a small portion of the 



phase space, "on the fly" ab initio molecular dynamics 
simulations can still be useful |Q . 

But even if an analytical representation of the ab ini- 
tio PES is available, there is still some computational 
effort to determine dissociation probabilities. It is a wide- 
spread believe that classical dynamical methods are much 
less time-consuming than quantum ones. This is cer- 
tainly true if one compares the computational cost of 
one trajectory to a quantum calculation. The quantum 
calculations, however, take advantage of the delocalized 
nature of the wave functions. A beam of particles ap- 
proaching the surface is described by a plane wave in the 
gas phase that includes all impact points. And a non- 
rotating molecule is represented by an j = rotational 
state, where j is the rotational quantum number and 
which contains all orientations with equal probability. 
Hence the averaging over initial conditions is done au- 
tomatically in quantum mechanics. Consequently, quan- 
tum dynamical simulations do not necessarily have to be 
more time-consuming that classical calculations, in par- 
ticular if one considers the fact that in wave-packet cal- 
culations the dissociation probability for a whole range 
of energies is calculated in one run or that in a time- 
independent coupled-channel method the sticking and 
scattering probabilites of all open channels are deter- 
mined simultaneously. 

The crucial difference between quantum and classical 
dynamics is that in the quantum dynamics the averaging 
is done coherently while it is done incoherently in the 
classical dynamics. For heavier molecules this is appar- 
ently not important, but for hydrogen dissociation this 
has important consequences, as will be shown below. 



III. HYDROGEN DISSOCIATION ON METAL 
SURFACES 

After introducing the theoretical concepts necessary to 
perform ab initio dynamics calculations, I will now dis- 
cuss applications. I start with the hydrogen dissociation 
on metal surfaces. 

While the dissociation of hydrogen at simple or noble 
metals is usually hindered by a substantial energy barrier, 
transition metal surfaces are rather reactive, i.e., already 
at low kinetic energies the dissociative adsorption prob- 
ability is often larger than 10 per cent. This reactivity 
has been attributed to the local density of states (LDOS) 
of the partly filled d-band at the Fermi level [119, 120 1 or 
close to the Fermi level [121,1221 of the transition met- 



als. Recent investigations have shown that the reactivity 
cannot be solely understood in terms of the LDOS at 
the Fermi level, but rather by the hybridization between 
molecular and metal states which involves the whole ri- 



band J45| jl23|jl24| ]. 

In noble metals the d-band is filled, which leads to 
the existence of a significant barrier towards dissociative 
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FIG. 1. Contour plot of the PES along a two-dimensional 
cut through the six-dimensional coordinate space of 
H 2 /Cu(lll) determined by DFT-GGA calculations. The in- 
set illustrates the orientation of the molecular axis and the 
lateral H2 center-of-mass coordinates, i.e. the coordinates 
X, Y, 6, and (j>. The coordinates in the figure are the H2 
center-of-mass distance from the surface Z and the H-H in- 
teratomic distance d. Energies are in eV per H2 molecule. 
The contour spacing is 0.1 eV. This cut corresponds to the 
minimum energy pathway (from ref. [^9|). 



adsorption of hydrogen. First I will discuss the ab ini- 
tio dynamics calculations for such an activated system, 
namely H2/CU which has been the benchmark system 
for the study of dissociative adsorption. Then I will ad- 
dress the dissociation on transition metals surfaces. For 
these latter systems ab initio dynamics calculations have 
proven their power by demonstrating the hitherto under- 
estimated efficiency of one particular dynamical mecha- 
nism, the steering. This chaper ends with a discussion 
of the influence of preadsorbates on the energetics and 
dynamics of hydrogen dissociation on metal surfaces. 



A. Dissociation on simple or noble metals 

The model system for the study of dissociative ad- 
sorption over the last years has been H2/CU. Many 
detailed experimental |125 12S ] as well as theoretical 
|,|J^j^,g,|^,|^,|l]j32|rdynamical investigations about 
this system exist. The early dynamical calculations were 
mostly performed on model potentials. Sometimes the 
used potential energy surfaces were derived from a small 
cluster calculation where the Cu surface had been mod- 



X[A] 

FIG. 2. The dissociation barrier height for H2/Cu(lll) as 
a function of the two lateral H2 center-of-mass coordinates 
X and Y determined by DFT-GGA calculations. The po- 
sitions of the copper substrate atoms are marked by "Cu". 
The molecular axis is parallel to the surface,its azimuthal 
orientation is illustrated by "H-H" (from ref. [39]). 



Then in 1994 there appeared two papers back-to-back 
in Physical Review Letters |39],flO) which represented the 
first detailed determination of a multi-dimensional PES 
as a function of all six hydrogen degrees of freedom by 
density- functional calculations. First of all these cal- 
culations showed that one has to go beyond the local 
density approximation (LDA) for the treatment of the 
exchange-correlation effects in density-functional theory 
calculations. LDA calculations yield an almost van- 
ishing barrier to dissociation although the experiments 
show that the height of this barrier should be about 
By including gradient corrections in 



0.5 eV ]125Hl29l. 

the exchange-correlation functional (the so-called gener- 
alized gradient approximation (GGA) Q) the general 
agreement with the experiment is greatly improved. Fig- 
ure [l] shows the GGA minimum barrier for H2/Cu(lll) in 
a so-called elbow plot. The PES is plotted as a function 
of the H2 center-of-mass distance from the surface and 
the interatomic H-H distance. The molecular orientation 
and the center-of-mass coordinates parallel to the surface 
are kept fixed. The calculations to determine the PES are 
not trivial; care has to be taken that the results are well- 
converged ^9,130|, and there are still some uncertainties 



elled by a Cu2 dimer |12C], or the PES was obtained from 
approximate methods like the effective medium theory 

BIB 



about the exact barrier heights [40 130 131]. GGA repre- 
sents a great improvement over LDA for certain systems 
as far as barrier heights are concerned, but it has its lim- 
itations |^5|,^], and the development of more accurate 
functionals is still going on (see, e.g., Ref. |86,132|). 
The two papers |3j||l| represented indeed a break- 



through in the determination of accurate and detailed po- 
tential energy surfaces. They challenged the gas-surface 
theorists with the fact that the interaction of hydrogen 
molecules with metal surfaces was much more corrugated 
than was anticipated before. Since the electronic density 



in front of metal surfaces is rather smeared out [ 1 19 , 133 1 , 
it was assumed that the barrier for dissociation depends 
only very weakly on the location within the surface unit 
cell. This picture of a fiat, structureless surface was seem- 
ingly confirmed by the experimentally found normal en- 



ergy scaling of the sticking probability [ 125 , 129 , 134 1 , i.e., 



the sticking probability is a function of the normal com- 
ponent of the incident kinetic energy alone. 

The ab initio calculations, however, showed that the 
barrier to dissociative adsorption for the molecular axis 
parallel to the surface varies by 0.7 eV within the surface 
unit cell (see Fig. |) |§. The bond-breaking process 
even on metal surfaces is a very localized process close to 
the surface. It involves the hybridization of molecular or- 
bitals with particular metal states, in particular d-states. 
These have a spatially strongly varying distribution re- 
flecting their symmetry properties. For an illustration 
of the local character of this interaction see for example 
Refs. [122 135 1 . Hence it is the chemical nature of the 
dissociation process that leads to its spatial variation. 

Motivated by the results of the ab initio calculations 
for the H2/CU PES, Darling and Holloway showed by 
using a three-dimensional model PES that strong corru- 
gation can be reconciled with normal energy scaling [ p2[ . 
For this to happen, the higher barriers have to be further 
away from the surface than the lower ones. This "bal- 
anced" corrugation ]3§| ] was indeed found in the ab initio 
calculations J3l|fl0fl , 

Shortly after the publication of the ab initio 
H2/Cu(lll) PES the first high-dimensional dynamical 
calculations were performed ||93|1 where the PES was en- 
tirely based on the first-principles calculations J3^]. A 
five-dimensional parametrisation of the ab initio PES 
has served as an input for the time-independent coupled- 
channel study. In this study the three center-of-mass co- 
ordinates and the interatomic spacing of the H2 molecule 
have been treated quantum mechanically, while the az- 
imuthal orientation of the molecule has been taken into 
account in a classical sudden approximation which works 
quite well in the H 2 / Cu system . 

These calculations demonstrated the importance of the 
multi-dimensionality for the determination of the stick- 
ing probability. Before, in the low-dimensional studies 
the barrier thickness has been used as a variable pa- 
rameter in order to reproduce the slope of the sticking 
curve. In low-dimensional calculations the width of the 
increase is basically caused by tunneling for energies be- 
low the barrier height and by quantum back scattering for 
energies above the barrier height. Furthermore, mixed 
quantum-classical calculations of the sticking probability 
for H 2 /Cu obtained by effective medium theory [^l|j32| 



showed only a modest dependence of the width of the 
sticking curve on the dimensionality of the calculations. 
The 5D calculations on the GGA-PES demonstrated that 
only in the tunneling regime the rise of the sticking proba- 
bility is determined by the barrier thickness J9^[ . Figure || 
shows the results for the calculated sticking probabilities 
of these 5D calculations. They are also compared to 2D 
calculations that included only the minimum energy path 
geometry plotted in Fig. [l| For energies above the min- 
imum barrier height the results of the 2D and the 5D 
calculations are very different. In this regime it is the 
whole barrier distribution that determines the width of 
the sticking curve. In the tunneling regime, on the other 
hand, the sticking curves of the 2D and the 5D calcula- 
tions are just shifted (Fig. ^p) demonstrating that in this 
regime it is indeed the minimum barrier thickness that 
governs the adsorption dynamics. 

As mentioned above, in the low-dimensional studies it 
was basically the barrier width that determined the rise 
of the sticking probability with kinetic energy. To re- 
produce the gradual increase of the experimental results, 
small barrier widths had to be assumed which were much 
smaller than the results obtain in the GGA calculations. 
In addition, to account for the experimental fact that vi- 
brational excitation of the impinging molecules is very 



efficient for increasing the sticking probability [ 125 1 , a 
late barrier after a strongly curved region according to 



the Polanyi's rules 13E] had been assumed in the model 
PES P]6|,|l0|] . The curvature of the reaction path in the 
GGA-PES (Fig. is much smaller than was previously 
assumed. It is not only the curvature, but also the low- 
ering of the vibrational frequency in the barrier region, 
which on the other hand was underestimated before, that 
caus es the vibrational effects in the dissociation process 

With regard to the effect of additional parallel momen- 
tum, Fig. y shows that for non-normal incidence the stick- 
ing probabilities fall upon the normal incidence sticking 
curve if they are plotted versus the normal kinetic energy 
for kinetic energies larger than 0.6 eV, i.e. in the regime 
where sticking is classically possible. The GGA-PES in- 
deed shows the topological features necessary for normal 
energy scaling (2^] , namely the balanced corrugation |3§| l 
with the larger barriers further away form the surface 
than the smaller. The 5D-calculations thus confirm the 
3D model calculations J2^] with respect to the role of the 
surface corrugation. 

However, for the tunneling regime Fig. ^b) reveals that 
additional parallel momentum enhances sticking. This is 
seemingly at variance with the molecular beam adsorp- 
tion experiments in which normal energy scaling has been 
found for all energies in the system H 2 / Cu( 1 11) [125,126]. 
However, this apparent discrepancy can be explained 
by the fact that at low kinetic energies sticking in the 
beam experiments is dominated by the vibrationally ex- 
cited molecules, and for these molecules the range where 
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FIG. 4. Classical and quantum dissociation probability as 
a function of kinetic energy for initially non-rotating H2 and 
D2 molecules under normal incidence at Cu(lll), determined 
by three-dimensional wave-packet calculations [Q. In these 
calculations the center-of-mass distance from the surface, one 
center-of-mass coordinate parallel to the surface and the polar 
orientation of the molecule were taken into account. 



potential. Due to the corrugation parallel momentum can 
be efficiently transferred to normal momentum for tun- 
neling particles, and this causes the enhancement of the 
sticking probability for additional parallel momentum in 
the low-energy regime. 

The GGA-PES for H 2 /Cu(lll) of Hammer and 
coworkers served also as an input for three-dimensional 
wave-packet calculations 94 in which, besides the H2 
center-of-mass distance from the surface, one H2 center- 
of-mass coordinate parallel to the surface and the po- 
lar orientation of the molecule were considered. The 
change of the vibrational frequency was taken into ac- 
count by keeping the molecule adiabatically in its vibra- 
tional ground state. This study particularly focused on 
the difference between the results of quantum and classi- 
cal calculations. Figure || shows the dissociation proba- 
bility as a function of the kinetic energy for energies close 
to the minimum barrier height. Classical and quantum 
results are shown for the two hydrogen isotopes H2 and 
D2. The difference between the two isotopes is due to 
the decrease in the zero-point vibrational energy during 
the dissociation process. Since the vibrational frequency 
of D2 is smaller by a factor of 1/V2 than the H2 fre- 
quency, for D2 less vibrational energy is transferred to 
the translational degree of freedom leading to the smaller 
dissociation probability. 

Rather surprising is the fact that for energies close to 
the dissociation barrier E ss 0.65 eV the classical prob- 
abilities are larger than the quantum ones. Naively one 
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FIG. 5. Sticking probability for H2 dissociation on Cu(100) 
determined by six-dimensional time-dependent quantum cal- 
culations on an ab initio PES (f rom Ref. ]5^|). The experi- 

The dotted line shows 




mental data are taken from Ref. [134 



the calculated probability for vibrational excitation in scat- 
tering. 



would expect that the quantum results should be larger 
due to tunneling which is absent in the classical calcula- 
tions. However, the PES is corrugated and anisotropic. 
For energies close to the minimum barrier energy only 
a limited range of molecular orientations and lateral co- 
ordinates are energetically accessible at the barrier posi- 
tion. The localisation of the wave function in the angu- 
lar and lateral degrees of freedom during the crossing of 
the minimum barrier region leads to the building up of 
zero-point energies in the quantum calculations. These 
zero-point energies cause an effective increase of the bar- 
rier height in the quantum dynamics which reduces the 
quantum dissociation probability compared to the classi- 
cal probability. The zero-point energies also cause quan- 
tization effects in the dissociation probability at higher 
energies which are visible as a weak step-like structure in 
Fig. |. 

Just recently the first six-dimensional wave-packet 
studies of the dissociation of H2 on Cu have appeared 
p2[ which confirmed the importance of the multi- 
dimensionality for the dissociation process. Fig. ^| shows 
the results for the sticking probability of H2 on Cu(100) 
p2[ . The calculations were based on a analytical fit to the 
PES obtained by density functional theory using GGA. 
The wave-packet method used a symmetry-adapted ba- 
sis set in the center-of-mass degrees of freedom parallel to 
the surface and for the molecular orientation. This causes 
a significant computational saving for initial conditions 
corresponding to normal incidence. 6D wave-packet cal- 
culations for non-normal incidence have not been per- 



formed yet. The theoretical results show some structure, 
in particular close to the threshold at 0.4 eV. These oscil- 
lations are attributed to resonances caused by the weak- 
ening of the molecular bond close to the surface [Q . 

The 6D sticking curve seems to be shifted by roughly 
0.2 eV to higher energies with respect to the experimen- 
tal results. A similiar shift is found in the 5D calculations 
for H 2 /Cu(lll) in Fig. §. This suggests that the GGA 
calculations overestimate the dissociation barriers. For 
H2/Cu(lll), convergence tests of the ab initio calcula- 
tions suggested that the barriers should be lowered by 
0.2 eV |3(|; the inclusion of the polar orientation, which 
was not considered in the 5D calculations, would shift 
the sticking curves to higher energies again, though. 

On the other hand, the sticking probability is directly 
measured only for kinetic energies up to approximately 
0.5 eV. H2 beams with higher kinetic energies can not be 
prepared in nozzle experiments. The experimental stick- 
ing probabilty for higher kinetic energies is derived from 
desorption experiments invoking the principle of detailed 
balance or microscopic reversibility. Because of this in- 
direct procedure, for example the saturation value of the 
sticking probability, i.e. the maximum value at high en- 
ergies, is not well-established experimentally; it has an 
uncertainty of a factor of two [126|. Hence it is not en- 
tirely clear yet, whether the discrepancy between theory 
and experiment is caused by experimental uncertainties 
or by approximations in the calculations. 

Another six-dimensional time-dependent dynamical 
study is devoted to the dissociation of H2 on Cu(lll) 
[ 139 1 . However, this study is not fully based on an ab 
initio PES. The authors claim that the available ab initio 
information is not complete enough, however, they have 
not collaborated with any group that is able to perform 
ab initio total-energy calculations. They use a LEPS po- 
tential for the parametrization of the H2 interaction with 
Cu(lll). The LEPS parametrization is considered by 
other authors jlTj to be not flexible enough to accurately 
describe the PES. 

Although already many detailed studies exist about 
the system H2/CU, experimentally as well as theoret- 
ically, there is still room for further explorations in 
this system. For example, current experimental stud- 
ies have focused on the rotational alignment of desorbing 
molecules | 140 - 142|| . These investigations address the is- 
sue of the anisotropy of the dissociation barriers. An- 
other noble metal system under current interest is the 
dissociation of hydrogen on silver. Detailed experimen- 



tal data exist for this system [143 144 1, and just recently 
also the PES for H2/Ag(100) was determined by first- 
principles calculations [135 1 . Hence there are still inter- 
esting results to be anticipated for the interaction of hy- 
drogen with noble or simple metal surfaces. 
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B. Dissociation on transition metals 



Molecular beam experiments of the dissociative ad- 
sorption of H2 on Pd(100) have found that the sticking 
probability initially decreases with increasing kinetic en- 
ergy of the beam | 145| , 146|| . Similiar results have been 
obtained for the H2 dissociation on various other tran- 
sition metal surfaces [145 147 152]. Such a behaviour 
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is usually found for atomic or non-dissociative molecu- 
lar adsorption in which the particles have to dissipate 
their excess energy to the substrate in order to remain 
on the surface pi. Since the energy transfer to the 
substrate becomes less efficient at higher energies, the 
sticking probability for atomic or molecular adsorption 
decreases with increasing energy. Consequently a sim- 
iliar process has been invoked to explain a decreasing 
sticking coefficient in dissociative adsorption, namely the 
precursor mechanism: before dissociation the molecules 
are temporarily trapped in a molecular adsorption state, 
the so-called precursor state. This trapping probability 
which decreases with increasing energy is supposed to be 
the rate-determining step. 

However, for hydrogen adsorption the large mass mis- 
match between adsorbate and substrate should make 
the energy transfer process inefficient. For the system 
H 2 /W(100)-c(2x2)Cu, e.g., it was shown (l5l) that for 
a hydrogen molecule impinging on a metal substrate the 
energy transfer to substrate phonons is much too small 
to account for the high sticking probability at low kinetic 
energies. 

Furthermore, density-functional theory calculations 
within the GGA found that there exist non-activated as 
well as activated paths to dissociative adsorption in the 
system H 2 /Pd(100), but the calculations gave no indi- 
cation of any molecular adsorption well, i.e. any pre- 
cursor state |4l],[4§. Two elbow plots of the H 2 /Pd(100) 
PES determined in these calculations are shown in Fig. ||. 
Fig. ^a) demonstrates that the dissociative adsorption of 
H 2 on Pd(100) is non- activated, i.e. reaction pathways 
without any hindering barrier exist. For the conditions 
of Fig. ||b), where the molecule approaches the surface at 
the on-top site, a barrier of approximately 0.15 eV exists. 
There seems to be a molecular adsorption well in front 
of the on-top site, but detailed calculations have shown 
that this well does not correspond to a local minimum 
of the PES (H). It is rather a saddle point of the PES 
in the multi-dimensional configuration space because the 
molecule can still follow a purely attractive path to disso- 
ciative adsorption if its center-of-mass degrees of freedom 
are allowed to relax. 

According to this theoretical results, the precursor 
mechanism can not be operative in the H 2 /Pd(100) sys- 
tem. As an alternative explanation it had been suggested 
that a decreasing sticking coefficient could be caused 
by a steering effect in a direct non-activated adsorption 
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FIG. 6. Contour plots of the PES 
along two two-dimensional cuts through the six-dimensional 
coordinate space of H^/Pd (100), so-called elbow plots, deter- 
mined by GGA calculations The coordinates in the 
figure are the H2 center-of-mass distance from the surface Z 
and the H-H interatomic distance d. The dissociation pro- 
cess in the Zd plane is illustrated in the inset. The lateral 
H2 center-of-mass coordinates in the surface unit cell and the 
orientation of the molecular axis, i.e. the coordinates X, Y , 
0, and cf) are kept fixed for each 2D cut and depicted above 
the elbow plots. Energies are in eV per H2 molecule. The 
contour spacing in a) is 0.1 eV, while it is 0.05 eV in b). 



process by King twenty years ago [153], but there had 
been no theoretical confirmation whether this mechanism 
could be efficient enough. In classical stochastic molec- 
ular dynamics simulations it was demonstrated how the 
anisotropy and corrugation could lead to a focusing of 
the dissociating molecules to reactive sites, but the en- 
ergy dependence of this process was not analysed |p8| . 
Furthermore, in two-dimensional dynamical treatments 
of t he H 2 /Pd(100) system no steering effect was observed 

The first dynamical evidence that the co-existence of 
non-activated with activated pathways to dissociative ad- 
sorption can lead to a initially decreasing sticking proba- 
bility was found in three-dimensional quantum dynamical 
calculations E6| . In this study a model PES with features 
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FIG. 7. Sticking probability versus kinetic energy for a 
hydrogen beam under normal incidence on a Pd(100) sur- 
face. Theory: six-dimensional quantum calculations for H2 
molecules initially in the rotational and vibrational ground 
state (dashed line) and with an initial rotational and en- 
ergy distribution adequate for molecular beam experiments 
(solid line) JS(J. H2 molecular beam ad sorp tion experiment 
under normal incidence (Rendulic et al. 145|): circles; H2 ef- 
fusive beam scattering experiment with an incident angle of 
of 6i = 15° (Rettner and Auerbach |146|): long-dashed line. 



derived from the GGA calculations of H 2 /Pd(100) @ 
was used. However, still large quantitative discrepancies 
to the experiment remained, as far as the relevant energy 
range in which steering is operative was concerned. 

Using a parametrization of the ab initio PES [ ^l|j42| 
based on 250 different configurations, six-dimensional dy- 
namical calculations of the dissociative adsorption and 
associative desorption of H2/Pd(100) have been per- 
formed |Q. These calculations were indeed the first in 
which all six degrees of freedom of the hydrogen molecule 
were treated quantum dynamically. Figure [?] presents the 
6D results of the sticking probability of H2/Pd(100) un- 
der normal incidence. These results are compared to the 
H2 molecular beam adsorption experiment by Rendulic, 
Anger and Winkler [ 145 1 . Furthermore, the experimen- 



tal results of Rettner and Auerbach for a reflection beam 
experiment under an angle of incidence of = 15° ]146| 
are plotted. 

First of all, a strong oscillatory structure is evident in 
the quantum results. Such oscillations have also been 
found in 4D wave-packet calculations of H2/W(100) |6l| l 
which is also a reactive system; those were, however, 
smaller. It has been shown that the energetic location of 
most of the peaks, in particular at low kinetic energies, 
can be related to the opening of new scattering channels 
with increasing kinetic energies, i.e. they correspond to 
threshold effects [154 , 155 1 . In the case of the diffrac- 
tion channels, the energetic position of the peaks is thus 



fully determined by the geometry of the surface. Their 
size depends sensitively on the symmetry of the initial 
conditions. For normal incidence they are much more 
pronounced than for non-normal incidence because the 
number of degenerate, symmetrically equivalent scatter- 
ing channels is reduced for non- normal incidenc e f31 ,155|. 
Holloway and coworkers have recently shown |15(| that 
the size of the oscillations also depends on the specific 
topology of the PES, in particular on the existence of 
possible resonance states in front of the surface. 

Rettner and Auerbach have searched for these oscil- 
lations hia^eflection experiment using a nearly effusive 
beam [ 146 157]. This kind of experiment allows a high 
energy resolution of the total scattering probability [ 158 1 . 
They have not found any indication of these oscillations. 
However, this is not an easy experimental task since sur- 
face impurities and the increasing hydrogen coverage dur- 
ing the scattering experiment destroy the coherence of 
the scattering event and suppress the oscillations. Fur- 
thermore, they have used an angle of incidence of 15° in- 
stead of normal incidence which also causes a suppression 
of the oscillations |^l| , since the symmetry of the initial 
conditions has a decisive influence on the magnitude of 
the oscillations. Simulated time-of-fiight (TOF) distribu- 
tions for H2 scattered at Pd(100) with the same angle of 
incidence as used in the experiment are very close to the 
experimental TOF distributions M] . The results of Ren- 



dulic, Anger and Winkler [14E] do not show any strong 
oscillatory structure either, but for an additional reason. 
An experimental molecular beam does not correspond to 
a monoenergetic beam in one specific quantum state. If 
one assumes an energy spread and a distribution of inter- 
nal molecular states typical for a beam experiment, the 
oscillations are almost entirely smoothed out in the 6D 
quantum results (solid line in fig. 0). 

The general agreement between theory and the exper- 
iments is satisfactory. As noted above, the ab initio PES 
does not have any molecular adsorption well, i.e., no pre- 
cursor state. In addition, due to the fact that the sub- 
strate is kept fixed, no energy transfer to the surface is 
taken into account in the calculations. Still the initial 
decrease of the sticking probability found in the experi- 
ments is well-produced in the quantum dynamical calcu- 
lations. Such an initial decrease of the sticking probabil- 
ity has also been found in classical molecular dynamics 
calculations on the same PES [§||L5§ and in quantum 
and classical calculations of the sticking probability of 
H2/W(100), too []si|| . Hence the underlying microscopic 
mechanism that is responsible for the decrease can also 
be investigated by analysing classical trajectories. This 
analysis yields that the large sticking probability at low 
kinetic energies is caused by the steering effect that be- 
comes less efficient at higher kinetic energies which then 
leads to the decrease of the dissociation probability. 

Figure || illustrates this general mechanism by showing 
snapshots of two ab initio molecular dynamics runs. For 
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a) Ei = 0.01 eV 




t = fs 




t = 160 fs 




t = 140 fs 




t = 215 fs 



b) £, = 0.12eV 




t = fs 




t = 80 fs 




t = 60 fs 




t = 165 fs 



FIG. 8. Snapshots of ab initio molecular dynamics trajectories for H2 molecules impinging on Pd(100) in order to illustrate 
the steering effect. The molecules are initially not vibrating and not rotating. For both trajectories the initial conditions are 
the same, except for the kinetic energy which is 0.01 eV in a) and 0.12 eV in b). 



both trajectories initially non- vibrating and non-rotating 
molecules impinge on the surface under the same condi- 
tions except for the initial kinetic energy which is 0.01 eV 
in Fig. ||a) and 0.12 eV in Fig. ||b). Far away from the 
surface the molecular axis is oriented almost perpendic- 
ular to the surface. In such an orientation the molecule 
cannot dissociate, the interaction with the surface is re- 
pulsive. The PES is, however, strongly anisotropic, and 
there are forces acting on the molecule to orient it to a 
parallel configuration to the surface. If the molecule is 
slow enough, as in Fig. ||a), it can indeed complete the 
rotation into this favorable orientation before hitting the 
repulsive wall of the potential. Close to the surface, after 
160 fs, the molecule has turned parallel to the surface and 
directly dissociates from this configuration. This process 
becomes less efficient at higher kinetic energies which is 
demonstrated in Fig. |8^>) where the initial kinetic energy 
is twelve times larger. Of course the same forces act upon 
the molecule, but now the molecule is too fast to be fully 
re-oriented. It hits the surface in an unfavorable con- 
figuration in which the interaction is repulsive. At the 
classical turning point there is a very rapid rotation cor- 
responding to a flip-flop motion, and then the molecule 
is scattered back into the gas-phase rotationally excited. 
By further increasing the kinetic energy, the molecules 
will eventually have enough energy to directly traverse 
the barrier region which causes the final increase in the 
sticking probability (see Fig. fjj). 

The events that are depicted in Fig. |[ direct disso- 
ciation or direct reflection, only correspond to the two 



limits of the possible outcome of a scattering event. A 
slow molecule approaching the surface can start rotat- 
ing caused by the anisotropy of the PES without directly 
dissociating. Due to the conversion of translational en- 
ergy into rotational energy it might lose so much trans- 
lational energy that it cannot escape into the gas phase 
again. The corrugation of the surface can also cause mo- 
tion parallel to the surface. In such a state the molecule is 
trapped into a dynamical precursor, which is not caused 
by energy transfer to the substrate but by energy transfer 
to rotational and kinetic energy of the molecule parallel 
to the surface. In ab initio molecular dynamics calcu- 
lations adsorption events have been found in which the 
molecule spent more than 5 ps in front of the surface be- 



fore dissociating [159]. Note that such events correspond 
to the existence of metastable states in quantum scatter- 
ing. These are very sensitive to the specific topology of 
the PES [l56l . 

Both the trapping into a molecular adsorption state 
by energy dissipation to the substrate and the suppres- 
sion of the steering effect lead to a sticking probability 
that decreases with increasing kinetic energy. It is there- 
fore desirable to find experimental properties that allow 
an unambiguous distinction between these two different 
mechanisms. As such a property the dependence of the 
sticking probability on the init ial r otational motion of the 
molecules has been suggested [ 16C ] . One of the properties 
of the PES that leads to steering, the anisotropy of the 
potential, also causes a suppression of the sticking proba- 
bility when the molecules are rapidly rotating. Molecules 
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orientationally averaged 
m i = 0, cartwheel rotation 
m, = j i: helicopter rotation 




1 2 3 

Initial rotational quantum number \. t 

FIG. 9. Sticking probability versus initial rotational quan- 
tum state ji for the system H2/Pd(100). Diamonds: ori- 
entationally averaged sticking probability (eq. |5^), triangles: 
rrii — (cartwheel rotation), circles: rrii = ji (helicopter ro- 
tation). The initial kinetic energy is Ei = 10 meV. 



with a high angular momentum will rotate out of a fa- 
vorable orientation towards dissociative adsorption dur- 
ing the time it takes to break the molecular bond. This 
is shown in Fig. |^ where the 6D sticking probability for a 
kinetic energy of E{ = 10 meV as a function of the initial 
rotational state is plotted. The diamonds correspond to 
the orientationally averaged sticking probability 
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The strong decrease with increasing initial rotational 
quantum number ji is caused by a suppression of the 
steering effect. Fast rotating molecules cannot be focused 
into a favorable orientation towards adsorption, they will 
be reflected from the surface. 

It is not easy experimentally to prepare a molecular 
beam in a single quantum state. Still, this rotational 
hindering of the steering effect has actually been con- 
firmed for H 2 /Pd(lll) |L61| , |l62j and also for H 2 /Pt(110) 
[ 163 1 . By seeding techniques the translational energy of a 
H2 beam can be changed in a nozzle experiment without 
altering the rotational population of the beam [161,162]. 
Also the different sticking probabilities of para- and n- 
hydrogen beams have been used to extract information 
about rotationally resolved results [163]. Figure [Tc| shows 
the results of a seeded beam experiment for H2/Pd(lll) 
by Beutl and coworkers [161]. The seeded beams have 
larger rotational energies than the unseeded beam at the 
same kinetic energies. It is said, they are rotationally 
hot. And indeed, especially at low kinetic energies the 
predicted strong suppression of the sticking probability 
by the rotational hindering is found, see Fig 
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FIG. 10. Experimental sticking probability versus kinetic 
energy for a H2 beam under normal incidence on a Pd(lll) 
surface. The seeded beams have larger rotational energies 
at id entical kinetic energies than the unseeded beam (after 



According to the principle of detailed balance, the 



suppresion of the sticking probability by the rotational 
hindering should be reflected by a population of rota- 
tional excited states in desorption which is lower than 
expected for molecules in thermal equilibrium with the 
surface temperature. This so-called rotational cooling 
has indeed been found for H2 molecules desorbing from 
Pd(100) [ 164, 165| and is also well reproduced by the six- 
dimensional quantum dynamical calculations |k| . 

In the precursor mechanism the molecules are assumed 
to be first trapped in a physisorption state. Such a 
state is caused by van-der-Waals forces which depend 
only weakly on the orientation of the molecule. Hence 
there are only small directional forces, and molecules ad- 
sorbed in such a physisorption state are able to rotate al- 
most freely Ji| [l66| . The trapping probability into the ph- 
ysisorption state and thus the sticking probability in the 
precursor model should therefore be almost independent 
of the initial rotational state, in contrast to the steering 
mechanism. In another series of experiments Beutl and 
coworkers have measured the dependence of the sticking 
probabilities on the initial rotational energies in the "clas- 
sical" precursor systems N 2 /W(100), CO/FeSi(100) and 
02/Ni(lll) [ 167 1 . They found no discernible influence of 
the rotational motion on the sticking coefficient in these 
systems. Hence the rotational dependence might indeed 
serve as a means to distinguish between the precursor 
mechanism and the steering effect. 

So far all the results presented were averaged over 
the azimuthal quantum number m, i.e, the results corre- 
sponded to orientationally averaged properties. Now the 
H2/Pd(100) PES is strongly anisotropic with regard to 
the molecular orientation. The most favourable configu- 
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ration towards dissociative adsorption is with the molec- 
ular axis parallel to the surface. Molecules that hit the 
surface in an upright position cannot dissociate, they are 
reflected back into the gas-phase. It is true that quan- 
tum mechanics does not allow for non-rotating, oriented 
molecules in the gas-phase, however, rotating molecules 
can show a preferential orientation. Fig. ^ also shows the 
effect of molecular orientation on the sticking probabil- 
ity. Molecules with azimuthal quantum number m = j 
have their axis preferentially oriented parallel to the sur- 
face. These molecules rotating in the so-called helicopter 
fashion dissociate more easily than molecules rotating in 
the cartwheel fashion (m = 0) with their rotational axis 
preferentially parallel to the surface since the latter have 
a high probability hitting the surface in an upright ori- 
entation in which they cannot dissociate. 

Experimentally it is hard to align a molecular beam of 
hydrogen. Again one can study the time-reverse process, 
the associative desorption. By laser-induced fluorescence 
(LIF) it is possible to measure the rotational alignment 
parameter A^\j) |16S| ], which is given by 



4 2) (i) 



3J' - J 1 



(23) 



(2) 

Aq (j) corresponds to the quadrupole moment of the 
oricntational distribution and assumes values of —1 

(2) 

Aq (j) < 2. Molecules rotating preferentially in the 

( 2^ 

cartwheel fashion have an alignment parameter A (j) < 
0, for molecules rotating preferentially in the helicopter 
fashion A^ (j) > 0. 

Experiments have in fact found a preferential pop- 
ulation of the helicopter mode in the desorption of 
D 2 /Pd(100) U|. Figure shows a comparison of these 
experimental results for D 2 with the six-dimensional cal- 
culation of the rotational alignment of H2 desorbing from 
a Pd (100) surface at a surface temperature of T s = 690 K 
|L70f . Experimentally there is no big difference in the 
alignment factors as a function of the rotational quan- 
tum number j between H2 and D2 [171|, hence the com- 
parison is meaningful. The agreement for j < 6 is quite 
satisfactory. Indeed the molecules desorb preferentially 
with their molecular axis parallel to the surface thus re- 
flecting the anisotropy of the PES. Due to computational 
restrictions the rotational alignment could only be calcu- 
lated for j < 6. For j = 7 and j = 8 the experiments 
show a vanishing alignment within the error bars which is 
rather surprising. It was argued that for these high rota- 
tional states an energy transfer from the rotations to the 
reaction coordinate could suppress the effect of the po- 
tential anisotropy ]169| |. Certainly these results deserve 
further clarification. 

The dependence of adsorption and desorptio n on ki 
netic energy, molecular rotation and orientation 
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FIG. 11. Rotationally alignment of hydrogen molecules 
desorbing fro m a Pd(100) surface. Boxes: experimental re- 
sults for D2 [16S|. Circ les: 6-D quantum dynamical calcula- 
tions for H 2 (from p7c|). 



angle of incidence |51|] , and the rot atio nally elastic and in- 
elastic diffraction of H2/Pd(100) [ 155 1 have been studied 
so far by six-dimensional ab initio dynamics calculations 
on the same PES. The results of these calculations have 
been compared to a number of independent experiments, 
and they are at least in semi-quantitative agreement with 
all of these experiments. Recently six-dimensional ab ini- 
tio quantum dynamical calculations of H2/Pd(100) have 
been performed |6£,172] which were based on a much 
larger set of DFT-GGA calculations jl35[ . This allowed 
a more detailed analytical representation of the PES and 
significantly improved the agreement between theoretical 
and experimental sticking probabilitities |66 172]. This 



shows that that ab initio dynamics calculations are ca- 
pable of adequately describing the hydrogen dissociation 
on transition metal surfaces. In addition, the difference 
between six-dimensional ab initio quantum dynamics and 
classical dynamics has been addressed [159]. As in the 



molecular vibration [ 137 1 , ro- vibrational coupling 
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case of the 3D-calculations of the hydrogen dissociation 
on Cu(lll) (see Fig. |J) the crucial difference between 
classical and quantum dynamics is not due to tunnel- 
ing, but to zero-point effects. However, caused by the 
higher dimensionality of the 6D-calculations the number 
of zero-point vibrational modes perpendicular to the re- 
action path is larger which makes the difference between 
classical and quantum results much more pronounced. 
This shows that a proper inclusion of zero-point effects 
is required in the treatment of hydrogen dynamics. I 
will address this issue in more detail in the next section 
about the hydrogen dissociation on an adsorbate covered 
surface. 

It is also well-known that palladium is able to absorb 
large amounts of hydrogen in the bulk and can there- 
fore be used for hydrogen storage |173|. The process of 
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how hydrogen might enter the bulk has been the issue of 
ab initio quantum dynamics studies employing two- and 
three-dimensional wave-packet calculations by Olsen et 
al. [174,175]. In order to determine the effect of surface 



motion on direct subsurface absorption of H 2 /Pd(lll) 
they have calculated a GGA-PES depending on two hy- 
drogen and one palladium degree of freedom ]175[ . They 
found that indeed the surface degrees of freedom play 
an important role at low kinetic energies as surface re- 
laxation can lower the effective barrier towards direct 
subsurface absorption substantially. Still this barrier 
is too large to explain the expe rimentally observed di- 
rect subsurface absorption |l76fl , probably due to the 
restricted dimensionality of these calculations, so that 
higher-dimensional calculations are desirable. 



Dissociation on an adsorbate covered metal 
surface 
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The presence of an adsorbate on a surface can pro- 
foundly change the surface reactivity. An understanding 
of the underlying mechanisms and their consequences on 
the reaction rates is of decisive importance for, e.g., de- 
signing better catalysts. Traditionally an "trial and er- 
ror" approach was used to improve the activity of a cata- 
lyst by adding some substances. Only recently this prob- 
lem has been addressed theoretically by ab initio methods 
[177). 

On Pd(100) it is experimentally well-known that the 
presence of sulfur leads to a large reduction of the hy- 
drogen dissociation probability [145,178|, it "poisons" 
the dissociation. While at the clean surface the disso- 
ciation probability is about 60% for a kinetic energy of 
E\ = 0.05 eV, at the sulfur-covered surface it drops below 



1% at the same energy 145 1. 

DFT-GGA calculations have shown that hydrogen dis- 
sociation on sulfur-covered Pd(100) is still exothermic, 
however, the dissociation is hindered by the formation of 
energy barriers in the entrance channel of the potential 
energy surface (PES) |43|] . Figure 12 shows two elbow 
plots of an analytical fit to the GGA-PES of H2 at the 
(2x2) sulfur covered Pd(100) surface for fixed molecu- 
lar orientation and lateral center-of-mass coordinates of 
the molecule. The minimum energy path is shown in 
Fig. |l2|a. The molecule approaches the surface above 
the fourfold hollow site. This is the site where the H2 
molecule is furthest away from the sulfur atoms on the 
surface. An analysis of the electronic structure calcula- 
tions has revealed that the building up of the barrier at 
this position is not related to a direct interaction between 
sulfur and the hydrogen molecule. Instead, it is caused 
by the sulfur-mediated downshift of the Pd d bands at 
the surface which leads to a population of anti- 

bonding molecule-surface states pa] . 



FIG. 12. Contour plots of the PES along two 
two-dimensional cuts through the six-dimensional coordinate 
space of H2/S(2x2)/Pd (100). The insets show the orienta- 
tion of the molecular axis and the lateral H2 center-of-mass 
coordinates, i.e. the coordinates X, Y, 8, and </>. The coor- 
dinates in the figure are the H2 center-of-mass distance from 
the surface Z and the H-H interatomic distance r. Energies 
are in eV per H2 molecule. The contour spacing in Fig. M2a 
is 0.1 eV, while in Fig. [f^j it is 0.2 eV. Fig. [l2^ , corresponds 
to the minimum energy pathway. 



Closer to the sulfur atoms the PES becomes strongly 
repulsive. This is illustrated in Fig. [l2|b. While the dis- 
sociation path over the Pd on-top position on the clean 
surface is hindered by a barrier of height 0.15 eV G3], the 
adsorbed sulfur leads to an increase in this barrier height 
to 1.3 eV. Directly above the sulfur atoms the barrier 
towards dissociation even increases to values larger than 
2.5 eV [Q for molecules oriented parallel to the surface. 
Here direct interaction between the sulfur atom and the 
H2 molecule is responsible for the repulsion and the high 
barrier [ ftij . 

On the analytical representation of the PES of H2 at 
the (2x2) sulfur covered surface six-dimensional dynam- 
ical calculations have recently been performed 103]. To 
my knowledge, this represents the most corrugated sys- 
tem of dissociative adsorption studied dynamically so far. 
Figure [l3| compares the quantum and classical results for 
the sticking probability as a function of the kinet ic en ergy 
of the incident H2 beam with the experiment [145]. In 
addition, also the integrated barrier distribution Pb(E), 



x cos 6d6 dcj) dX dY (24) 
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FIG. 13. Sticking probability versus kinetic energy for a 
H2 beam under normal incidence on a S(2 x2)/Pd(100) sur- 
face. Full dots: experiment (from ref. | 145 |); Dashed-dotted 
line: Integrated barrier distribution, which corresponds to the 
sticking probability in the hole model jnj; Solid line: Quan- 
tum mechanical results for molecules initially in the rotational 
and vibrational ground-state; Dashed line: Classical results 
for initially non-rotating and non-vibrating molecules. The 
inset shows the quantum and classical results at low energies. 



is plotted. Here 9 and <f> are the polar and azimuthal 
orientation of the molecule, X and Y are the lateral co- 
ordinates of the hydrogen center-of-mass. A is the area 
of the surface unit cell. Each quadruple defines a cut 
through the six-dimensional space (see Fig. ^| for exam- 
ples) , and Eb is the minimum energy barrier along such a 
cut. The function is the Heavyside step function. The 
quantity Pb(E) is the fraction of the configuration space 
for which the barrier towards dissociation is less than E; 
it corresponds to the sticking probability in the classical 
sudden approximation or the so-called "hole model" jl7) . 

The calculated sticking probabilities are significantly 
larger than the experimental results. The onset of dis- 
sociative adsorption at E\ w 0.12 eV, however, is well 
reproduced by the calculations. This onset is indeed also 
in agreement with the experimentally measured mean ki- 
netic energy of hydrogen molecules desorbing from sul- 
fur covered Pd(100) [179|. The theoretical results are 
relatively insensitive to variations in the barrier heights 
of the order of 0.2 eV [103]. It might therefore well be 
that uncertainties in the experimental determination are 
responsible for the difference. The exact sulfur cover- 
age in the experiment was not very well characterized. 
The sulfur adlayer was obtained by simply heating up 
the sample which leads to segregation of bulk sulfur at 
the surface. The sulfur coverage was monitored through 
the ratio of the Auger peaks Si32/Pd33o [145]. The set 
of experimental data shown in Fig. [l^ were obtained 
for half of the maximum segregation coverage of sulfur 



mode 


ZPE (eV) 


H-H vibration 


0.253 


polar rotation 


0.016 


azimuthal rotation 


0.013 


translation perpendicular to molecular axis 


0.027 


translation parallel to molecular axis 


0.027 


sum 


0.336 



TABLE I. Zero-point energies (ZPE) of the H2 molecule at 
the minimum barrier position. The H2 configuration corre- 
sponds to the situation of Fig. |l2]a. The gas-phase zero-point 



energy of H2 is ~Hu} gas = 0.258 eV. 



on Pd(100). This saturation coverage is assumed to be 
roughly 65 rs 0.5. Apart from uncertainties in this de- 
termination of the sulfur coverage also subsurface sulfur 
might have influenced the measurements. In ref. [178] 
a linear decrease of the hydrogen saturation coverage 
with increasing sulfur coverages was found. At a sul- 
fur coverage of 6s = 0.28 hydrogen adsorption should be 
completely suppressed. This infers a site-blocking effect 
which is at variance with the ab initio calculations of the 
PES for this system 



43 44 1 . However, these seemingly 
contradicting results and also the discrepancy between 
calculated and measured sticking probabilities could be 
reconciled if subsurface sulfur (which is not considered in 
the calculations but which might well be present in the 
experimental sample) has a decisive influence on the hy- 
drogen adsorption energies. It is certainly desirable that 
the effect of subsurface sulfur on the hydrogen adsorption 
in this system will be investigated. 

Figure |l3| shows furthermore that the classical molecu- 
lar dynamics calculations over-estimate the sticking prob- 
ability of H 2 at S(2x2)/Pd(100) compared to the quan- 
tum results. There are two important quantum phe- 
nomena not taken into account by classical calculations: 
tunnelling and zero-point effects. For energies smaller 
than the minimum barrier height E^ the classical sticking 
probability is of course zero, whereas the quantum results 
still show some dissociation due to tunneling, as the inset 
of Fig. [l^ reveals. But for energies larger than E\, the 
classical sticking probability rises to values which are al- 
most 50% larger than the quantum sticking probabilities. 
Tunnelling increases the quantum transmission probabil- 
ity with respect to the classical result, hence tunnelling 
would have the opposite effect as observed in Fig. [l^. 
Thus only zero-point effects can be responsible for the 
large difference [[31,159 



Since the hydrogen molecular bond is still almost in- 
tact at the minimum barrier position, the H-H vibra- 
tional frequency should not be too different from its gas- 
phase value there. However, the minimum barrier is lo- 
calized in the lateral coordinates parallel to the surface 
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FIG. 14. Illustration of steering at the (2x2) sulfur covered 
Pd(100) surface. Center-of-mass trajectories of H2 molecules 
are shown impinging with an kinetic energy of Ei = 0.3 eV on 
the diagonal of the (2x2) unit cell. Dashed lines correspond to 
reflection events while the full lines depict adsorption events. 



FIG. 15. Dependence of the quantum sticking probability 
versus ikinetic energy for a hydrogen beam under normal inci- 
dence at a S(2x2)/Pd(100) surface on the initial vibrational 
state v. Solid line: Vi = 0, dash-dotted line: m = 1. The 
molecules are initially in the rotational ground state. 



and in the rotational degrees of freedom of the molecule. 
The associated modes become "frustrated" which leads 
to zero-point energies in these degrees of freedom. This 
is confirmed by Tab. | where the zero-point energies of 
the hydrogen molecule at the minimum barrier position 
are collected. Note that the zero-point energies in the 
frustrated rotational and lateral modes are still small, in 
particular the rotational modes. This is caused by the 
fact that the molecular bond is essentially not elongated 
at the minimum barrier position so that the molecular 
interaction with the surface is still rather isotropic. 

However, the sum of all zero-point energies at the min- 
imum barrier position is in fact 0.08 eV larger than the 
gas-phase zero-point energy of hydrogen which is just half 
the molecular vibrational frequency ^huj gas — 0.258 eV. 
The combined effect of all zero-point energies at the min- 
imum barrier position is to enlarge the effective barrier 
height by 0.08 eV to E e b ff = 0.17 eV. The absence of 
zero-point effects in the classical dynamics causes the 
strongly enhanced sticking probability in the classical cal- 
culations. Still there is already a significant quantum 
sticking probability for kinetic energies below this effec- 
tive barrier height. A detailed dynamical analysis sug- 
gests that this is d ue to the combined effect of steering 
and tunneling [ 180 1 . Classically steering does not help to 
traverse the barrier region for energies below the barrier 
height. Also for energies only slightly larger than the bar- 
rier height steering is not very efficient. This is different 
in quantum dynamics where steering can already occur 
in the tunneling regime. Therefore it is very effective for 
energies close to the minimum barrier height with respect 
to the classical results. 



Although the dissociation of hydrogen on the sulfur 
covered Pd(100) surface in contrast to the clean surface 
is hindered by barriers, there is also significant steering 
of the impinging molecules to low-barrier configurations 
operative. This is reflected by the fact that the calcu- 
lated sticking probabilities are much larger than what one 
would expect from the hole model [|l7| (see Fig. |l3|). This 
steering has been confirmed by analysing swarms of tra- 
jectories and is illustrated in Fig. |l4|. The center-of-mass 
trajectories of H2 molecules impinging on the diagonal 
of the unit cell of the (2x2) sulfur covered Pd(100) sur- 
face are shown. The initial kinetic energy is 0.3 eV. The 
molecular axis is parallel to the surface, but the initial 
azimuthal orientation is chosen in such a way that in the 
sudden approximation all initial conditions do not lead 
to dissociative adsorption since all incoming molecules 
hit barriers larger than 0.3 eV. The strongly repulsive 
potential directly above the sulfur atoms extends rather 
far into the gas-phase [Q. Incoming molecules are thus 
diverted from this high-barrier sites. Near to the sul- 
fur atoms the molecules are still scattered back into the 
gas-phase, but closer to the center of the (2x2) unit cell 
the molecules are focused to the minimum barrier site. 
This focussing is accompanied by a reorientation of the 
molecular axis so that the molecule dissociates in the 
configuration depicted by the inset of Fig. |l2|a. 

The effect of initial vibrational motion on the disso- 
ciation process is demonstrated in Fig. [l5| where state- 
specific sticking probabilities as a function of the inci- 
dent kinetic energy are plotted. As Fig. [l~2| a shows, the 
minimum barrier to dissociation is at a position where 
the bond-length of the molecule is still not significantly 
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elongated. In such a situation it is usually anticipated 
that the vibrational and translational degrees of free- 
dom are almost uncoupled so that vibrational energy 
of the impinging molecules cannot be used to overcome 
the barrier (see, e.g., Ref. Q). Consequently, it has 
been predicted |43) that the sticking probability of H 2 at 
S(2x2)/Pd(100) should show no strong dependence on 
the initial vibrational state of the molecule. This pre- 
diction corresponds to the so-called Polanyi rules which 
have been formulated for gas-phase dynamics thirty years 
ago |136[ . Apparently the rules have to be modified for 
strongly corrugated surfaces. Figure |l^ demonstrates 
that initial vibrational excitation leads to a significant 
increase in the sticking probability in this early barrier 
system. As Fig. |l2| b reveals, there are molecular configu- 
rations for which the bond is significantly extended at the 
barrier position. Adsorbing molecules do not only probe 
the minimum barrier of the PES, but also other regions 
of the PES where vibrational energy can be efficiently 
used to overcome the dissociation barriers. This result 
of the vibrationally enhanced dissociation is actually in 
agreement with the experimentally observed vibrational 
over-population in thermal hydrogen desorption from sul- 
fur covered Pd(100) 1 181 1 invoking the principle of micro- 
scopic reversibility. 

It is interesting to note that the quantum sticking 
probability for molecules in the first excited vibrational 
state is close to the classical sticking probability for non- 
vibrating molecules (see. Fig. |l^). For the hydrogen 
dissociation on the clean Pd(100) surface is has been 
shown that zero-point effects can cancel out of the dy- 
namics calculations if the sum of all zero-point energies 
remains constant along the dissociation path ]159[ . For 
H2 at S(2x2)/Pd(100) this sum increases as Tab. D] shows. 
Now the H-H vibration corresponds to the fastest mode 
in the dissociation dynamics, and a comparison of 5D vi- 
brationally adiabatic with full 6D calculations has shown 
that the vibrations indeed follow the change of the vibra- 
tional frequency during the dissociation almost adiabati- 
cally 1 137]. For molecules in the first excited vibrational 
state the transfer from vibrational energy to translational 
energy due to the lowering of the vibrational frequences 
can compensate the increase in the zero-point energies in 
the rotational and lateral modes leading to the described 
cancellation effect. 



IV. HYDROGEN INTERACTION WITH 
SEMICONDUCTOR SURFACES 

The ab initio studies of hydrogen dissociation on metal 
surfaces provide a rather complete picture of the disso- 
ciation process. The quantitative agreement can still be 
improved by, e.g., the inclusion of dissipation processes 
to substrate phonons or electron-hole pairs, the general 
mechanisms, however, seem to be understood (of course 



one should always be careful with such optimistic state- 
ments). For the interaction of hydrogen with semicon- 
ductor surfaces, on the contrary, there is still no general 
agreement on the fundamental mechanism, especially in 
the well studied system H2/SL The reason lies in the 
fact that the rearrangement of the semiconductor surface 
upon hydrogen adsorption can no longer be neglected as 
for metal systems. This makes the adsorption/desorption 
dynamics and their theoretical description more compli- 
cated. 

In the following I will first discuss the dissociation dy- 
namics of hydrogen on Si(100). Then I will address a 
recent ab initio study of the reaction of atomic hydrogen 
with the hydrogen-passivated Si(100) surface. 



A. Dissociation on semiconductor surfaces 

The system H^/Si plays the same role for the hydrogen 
dissociation on semiconductors that H2/CU has for metal 
surfaces: it serves as a model system for which an abun- 
dance of experimental and theoretical studies exist (see 
for example the recent review by Kolasinski fl56|). Be- 
sides, this system is also of great technological relevance 
for, e.g., the etching and passivation of Si surfaces or the 
growth of Si crystals. 

Desorption of hydrogen from Si(100) shows first-order 
kinetics [182-185]. For associative desorption one nor- 
mally expects second-order kinetics since two atoms have 
to find each other on the surface before they can desorb. 
The unusual first-order desorption kinetics has been ex- 
plained by a prepairing mechanism ]91 , 183| ,184]: Desorb- 
ing molecules originate from the same dimer since it is 
energetically favorable for two hydrogen atoms to bind on 
the same dimer rather than on two independent dimers 
(see Fig. |l^a). 

The so-called barrier puzzle has further fueled the in- 
terest for this system: While the sticking coe fficient of 
molecular hydrogen on Si surfaces is very small 1 186 189 
indicating a high barrier to adsorption, in a laser-induced 
desorption experiment an almost t herm al mean kinetic 
energy of the molecules was found [ 190 1 indicating a low 
adsorption barrier. In order to explain this puzzle it was 
suggested to take the strong surface rearrangement of Si 
upon hydrogen adsorption into account [ 190 1 : The hy- 
drogen molecules impinging on the Si substrate from the 
gas phase typically encounter a Si configuration which 
is unfavorable for dissociation, while desorbing hydrogen 
molecules leave the surface from a rearranged Si configu- 
ration with a low barrier. In the adsorption process such 
a surface rearrangement can only be achieved by thermal 
excitations of the lattice since due to the mass mismatch 
the silicon substrate atoms are too inert to change their 
configuration during the time the hydrogen molecule in- 
teracts with the surface. The strong surface rearrange- 
ment of Si(100) upon hydrogen adsorption/desorption is 
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FIG. 16. a) Hydrogen covered Si (100) surface (monohydride) . b) Snapshots of a trajectory of D2 desorbing from Si (100) 
starting at the transition state with the Si atoms initial at rest pq| . The dark Si atoms correspond to the Si positions after the 
desorption event, c) Clean anti-buckled Si (100) surface. 



illustrated in Fig. [L6L 

With the assumption of a lattice rearrangement en- 
ergy of about 0.7 eV existing experimental adsorption 
and desorption results could be reproduced by quantum 
dynamical model calculations |9,191,192|. The predicted 
strong dependence of the adsorption probability on the 
surfa ce temperature was later confirmed experimentally 
[188J. 

So far the dissociation mechanism seemed to be at 
least qualitatively understood. However, there are some 
disturbing experimental uncertainties. Measurements 
of the dissociative adsorption probability of H2/Si(100) 
differ by almo st three orders of magnitude from each 
other 187 18£]. In order to e xplain the large difference 
it had been speculated [189], e.g., th at th e Si sample 
used for the adsorption study of Ref. [187] exhibits an 
unusual low barrier due to its high dopand concentration 
of n w 10 19 cm" 3 . 

From the theoretical point of view, the situation is 
even more confusing. The exact adsorption and des- 
orption mechanism is strongly debated. Total-energy 



surface was found to reconstruct with a symmetric dimer 



calculations using the cluster approach [193-198] have 
found activation barriers for associative desorption of 
hydrogen from Si(100) which are roughly 1 eV larger 
than the experimentally found value of about 2.5 eV 



[ 183, 184 1 . Based on these findings, the prepairing mech- 
anism has been disregarded, and defect-mediated des- 
orption mechanisms had been proposed [193, 195 - [200 1. 
These mechanisms, however, have difficulties explaining 
the observed first-order desorption kinetics and the mea- 
sured high prefactor of the desorption rate (for a thor- 
ough discussion of these calculations, see the review by 
Doren (9l]]). The question arises whether a cluster really 
represents a good model for an infinite substrate. For 
example, in most of the cluster calculations the Si(100) 



structure [195 199 201] while in slab calculations the anti- 
symmetric buckled dimer structure is the lowest energy 
configuration |202| , |203| which is also well-established ex- 
perimentally |204| " |206] ] . It was suggested that the re- 
laxation constraints used in the cluster studies do not 



model the Si(100) surface correctly [207]. Furthermore, 
the reported barriers for the same process differ by almost 
1 eV between three different sets of cluster calculations 



[193 -195 1 indicating the difficulty of such calculations. 

On the other hand, in detailed DFT-GGA calculations 
of the H2/Si(100) pote ntia l ene rgy surface (PES) using 
the superc ell a pproach | 208j -211| and also in clusters cal- 
culations |212| desorption barriers for hydrogen atoms 
originating from the same Si dimer were obtained that 
are in agreement with the experimentally determined val- 
ues. This pre-pairing mechanism, as mentioned above, is 
consistent with the observed first-order kinetics in ex- 



periment [182-185]. The adsorption barriers were found 
to be about 0.3 - 0.4 eV with the substrate rearrange- 
ment energy being roughly 0.15 eV. These barrier heights 
are also controversely discussed. It was argued that the 
GGA-PW91 functional used in the DFT studies under- 
estimates the H2 elimination barriers from silancs, i.e., 
small molecules pTj. It is, however, questionable whether 
the results for small molecules are directly transferable 
to surface systems. 

Subsequently, based on the information from the slab 



calculations [208, 209 two different three-dimensional 
quantum dynamical studies were performed [ 213 , 214 1 , 



where the relaxation of the Si substrate upon hydrogen 
adsorption was approximated by one idealized Si phonon 
coordinate. Although both calculations used the same 
ab initio energies as a source, the results of the dynam- 
ical calculations did not agree quantitatively. However, 
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FIG. 17. Ab initio molecular dynamics trajectory of D2 desorbing from Si(100). The initial conditions correspond to a surface 
temperature of T s — 920 K. a) Thermal desorption trajectory started close to the transition state (TS). Upper panel: x and 
z coordinates of the trajectories of the two deuterium atoms and the two Si atoms of the dimer underneath the TS. The final 
positions of the atoms are marked by the large dots. Lower panel: energy redistribution as a function of the last 140 fs of the 
run-time of the trajectory. The different energy channels are described in the legend of the figure, b) Time-reversed trajectory 
of the desorption event started at the TS. The upper and lower panels correspond to the ones in a). 



in both studies the dynamical coupling of the desorp- 
tion path to the Si substrate was low, i.e., the desorb- 
ing molecules had more excess translational energy than 



found in the experiment 19C]. 

These different uncertainties together with the ob- 
served unusual first-order desorption kinetics make the 
H2/Si(100) system to one of the most controversely dis- 
cussed systems in the field of gas-surface dynamics. Since 
barrier calculations alone will not settle the controversy 
about the desorption mechanism [fjl| , dynamics calcula- 
tions are necessary in order to find out how the potential 
energy at the barrier is distributed over the various de- 
grees of freedom of this system (hydrogen vibration, ro- 
tation, and translational energy, and vibrations of the Si 
substrate). Due to the importance of lattice relaxation 
effects on the desorption dynamics, substrate degrees of 
freedom have to be taken into account in the dynamics 
simulations. This makes a full quantum dynamical treat- 



ment not possible at the moment. Ab initio molecular 
dynamics calculations using both the supercell approach 
p8[ and clusters |^| have therefore been performed. In 
the supercell calculations DFT-GGA was employed to de- 
termine the desorption dynamics of deuterium desorbing 
from the monohydride (MH) state using a (2 x2) surface 
unit cell and a five-layer Si slab (see Fig. |ll|) . In the 
cluster calculations hydrogen desorption was treated on 
the complete active space self-consistent field (CASSCF) 
level. Both desorption from a single monohydride us- 
ing a SigHi2H2 cluster and from an isolated dihydride 
(ID) using a SiioHi4H 2 cluster were studied where 
H denotes the modified hydrogen atoms terminating the 
cluster. 

Since the barrier to associative desorption of hydro- 
gen from Si(100) is rather high, in both sets of ab initio 
molecular dynamics calculations [ j58|]59| the trajectories 
were not started with the hydrogen at the atomic adsorp- 
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tion positions. Instead, the molecular dynamics calcula- 
tions were initialized with the hydrogen atoms close to 
the transition state with initial velocity distributions that 
correspond to the surface temperatures used in the ex- 
periments. In the supercell calculations, additionally the 
time-reversed trajectories starting at the transition state 
were determined in order to make sure that all trajecto- 
ries correspond indeed to desorption and not to scattering 
events @. 

In the DFT-GGA calculations additional trajectories 
were determined with the Si lattice initially at rest, i.e., 
at a surface temperature of T s = K. Snapshots of such a 
calculated trajectory are shown in Fig. |l6|b). The dark Si 
atoms correspond to the relaxation of the Si lattice after 
the desorption event. Approximately 0.1 eV of the poten- 
tial energy at the transition state is transferred to vibra- 
tions of the Si lattice which is a rather large amount com- 
pared to hydrogen/metal systems. At the transition state 
the interatomic hydrogen distance is about 40% larger 
than the gas-phase bond length; consequently molecular 
vibrations are excited during the desorption, as is well 
known for a long time in associative desorption studies 
(see, e.g., ref. @). 

Similar results have been found in the DFT-GGA cal- 
culations for initial conditions corresponding to a surface 
temperature of T s = 920 K. In total 34 "thermal" des- 
orption trajectories were calculated. Figure |l7]a) shows 
an example of such a trajectory. First the surface was 
equilibrated for approximately 1 ps with the deuterium 
atoms kept close to the transition state by an auxiliary 
potential. Then the extra potential was switched off, 
the molecule was allowed to desorb and the distribution 
of the energy into the various degrees of freeedom was 
monitored. 

The projection of the trajectories of the desorbing deu- 
terium molecule and of the Si dimer closest to the transi- 
tion state onto the xz-plane is shown in Fig. |l7]a). Due to 
the asymmetry of the transition state the molecule does 
not desorb in the direction normal to the surface. The os- 
cillations in the trajectories correspond to the vibrational 
excitation of the desorbing D2. In addition, the energy 
redistribution during the last 140 fs of the run is plotted 
where, e.g., the vibrational excitation can be clearly fol- 
lowed. The quenching of the rotational motion is evident 
which is caused by the strong anisotropy of the PES. Fur- 
thermore, the deuterium molecule is strongly accelerated 
during the desorption event and ends up with a kinetic 
energy of roughly 0.4 eV which is larger than the value 
of IksTs = 0.16 eV expected for thermal equilibrium. 

The configuration at which the extra potential was 
switched off was also used to start a time-reversed tra- 
jectory. As mentioned above, the main reason was to 
check whether the trajectory describes a true reactive 
event. On the other hand, the time-reversed trajectory 
simultaneously corresponds to a dissociative adsorption 
event and allows to follow the energy distribution dur- 



ing the adsorption. Once the molecule has passed the 
transition state, the single atoms are strongly acceler- 
ated towards their atomic adsorption positions; together 
they can gain an energy of 2.5 eV. Figure |l7|b) shows 
that indeed the two D atoms first become very fast, in 
particular the left deuterium atom that has the longer 
distance to its adsorption position. Within 50 fs the D 
atoms have gained about 1.5 eV in kinetic energy. Al- 
though the mass mismatch between the deuterium and 
the silicon atoms is so large, there is still a significant en- 
ergy transfer to the silicon lattice due to the high kinetic 
energy of the D atoms. They bounce into the lattice so 
violently that within 150 fs 1 eV of the kinetic energy is 
transferred to the silicon atoms, as Fig. |l7|b) reveals. Fur- 
thermore, Fig. [ll] nicely demonstrates how the atoms of 
the Si dimer are moving towards a symmetric configura- 
tion in the case of the D2 dissociative adsorption, while in 
the case of the associative desorption the buckling angle 
is further increased with respect to the transition state 
configuration. 

The mean total, kinetic, vibrational, and rotational 
energies and final angles of the hydrogen molecules des- 
orbing from Si(100) obtained from the ab initio molecu- 
lar dynamics studies are compared to the experimental 
results in Tab. || Since there is apparently some confu- 
sion about doing a reasonable comparison between the- 
ory and experiment, I will briefly expound on the given 
values. The vibrational and rotational population of hy- 
drogen molecules desorbing from Si(100) was measured 
state-specifically using resonance-enhanced multiphonon 



ionization detection 215]. This gives information about 
the population of the quantum mechanical eigenstates of 
the hydrogen molecule. For example, it has been mea- 
sured that at a surface temperature of T 8 = 780 K 1% of 
desorbing H2 molecules are in the first vibrational eigen- 
state [ 215 1 . This corresponds to a fraction that is 20 
times larger than expected in thermal equilibrium and is 
described by the term "vibrational heating" . In Ref. |59fl 
this population was used to determine the mean vibra- 
tional energy of desorbing H2 molecules. They obtained 
a value of 0.264 eV. However, 98% of this value comes 
from the zero-point energy of the ground-state molecules, 
only 2% from molecules in the first excited vibrational 
state. Had the experiments not given any vibrational 
heating, this mean vibrational energy would have hardly 
been changed at all. In the ab initio molecular dynam- 
ics calculations the dynamics of the hydrogen nuclei is 
treated fully classically, there are no zero-point ener- 
gies considered. Since the H2 zero-point energy is so 
large (three times larger than the thermal energies in 
the experiment) it makes no sense comparing a value de- 
rived from a quantum mechanical analysis with a clas- 
sically obtained value. Therefore I have put in Tab. || 
for the experimental mean vibrational energy the value 
(E V ib) exp = ksTvii,, where T v n, = 1700 K was taken from 



the measured population 1 215] . This procedure is also not 
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(Etot) 


(Ekin) 


{E V ib} 


(E ro t) 


(Of) 


D 2 , MH, DFT-GGA 


0.70 ±0.18 


0.55 ±0.14 


0.12 ±0.09 


0.03 ±0.05 


39.6° ± 10.1° 


H 2 , MH, CASSCF 


0.72 ±0.20 


0.6 ±0.1 


0.09 ±0.08 


0.03 ±0.02 


56° ± 7° 


H 2 , ID, CASSCF 


0.94 ±0.20 


0.7 ±0.1 


0.2 ±0.1 


0.04 ±0.03 


43° ± 17° 


D 2 , Experiment 


0.34 ±0.05 


0.17 ±0.03 


0.15 ±0.03 


0.03 ±0.01 


28° - 31° 



TABLE II. Mean energy distribution and final desorption angle of hydrogen molecules desorbing from a Si(100) surface. 
The results of the DFT-GGA calculations are averaged over 34 desorption trajectories corresponding to D 2 desorption from 
the monohydride (MH) phase at a surface temperature of T a = 920 K (fcsT s = 0.079 eV) |Es| |. The cluster calculations at the 
CASSCF level for H 2 desorption from a MH and an isolated dihydride (ID) are averaged over ten trajectories for each case 
at a surface temperature of T a — 780 K (fcsT s = 0.067 eV) j5i|. The experimental results are collected from Ref. (56). The 
translational energy has bee n be en determined in a laser-induced desorption experiment where the surface temperature was 
estimated t o be T s — 920 K |190| while the vibrational and rotational energy have been measured at a surface temperature of 



T a — 780 K [215 . The mean vibrational energy corresponds to the classical value using the vibra tional temperature measured in 
the experiment (see text). The experimental mean desorption angles are derived from Ref. [ 216 1 , where the angular distribution 
in desorption has been measured for two different hydrogen coverages. 



ideal, but certainly closer to a classical analysis. A real 
quantitative comparison between experiment and theory 
can only be done for a quantum treatment of the hy- 
drogen vibrational dynamics. In the three-dimensional 
quantum dynamical studies, e.g, the correct experimen- 
tal population of vibrationally excited states in desorp- 
tion is reproduced [191 .192 , 213 1 . These problems do not 
arise for the mean translational and rotational energies. 
Since these modes are free in the gas-phase, there are no 
zero-point energies associated with them. 

Furthermore, in Ref. [ p^ a mean experimental kinetic 
energy of 0.34 eV for H2 desorption was derived which 
was said to be scaled from the experimental value of 
0.17 eV for D2. This scaling seems very questionable 
to me. H2 and D2 have exactly the same interaction po- 
tential. The gain in kinetic energy is determined by the 
change in the potential energy which is the same for H2 
and T>2- And classically different isotopes follow exactly 
the same trajectories as a function of the kinetic energy 
if no energy transfer to substrate degrees of freedom is 
considered Hence H2 and D2 should have approxi- 
mately the same kinetic energy in desorption. There can 
be an isotope effect due to quantum mechanical effects 
or due to energy transfer to the silicon lattice. However, 
the sign and magnitude of this isotope effect is not clear 
a priori; it has to be determined by detailed calculations. 
Certainly it does not follow such a simple scaling law as 
assumed in Ref. |5j| . 

Comparing the theoretical and experimental data in 
Tab. EL it should furthermore be noted that the reported 
theoretical values do not really correspond to any ther- 
mal average. For that purpose the number of 34 and 
10 calculated trajectories, respectively, is much too low. 
Hence the theoretical results have to be interpreted cau- 
tiously. Still they should be able to indicate the general 
trends in the energy distribution of desorbing molecules 
since in thermal desorption the molecules mainly origi- 



nate from a small fraction of the configuration space close 
to the transition state (the so-called keyhole effect |^6| ) . 
As for the mean desorption angle, note that this value 
correspond to an average over all azimuthal angles with 
< Of < 90°. Hence even if the angular distribution for 
fixed azimuth is peaked in the direction normal to the 
surface, (Of) will still be larger than zero. 

It is obvious that all simulations reproduce the vibra- 
tional heating, i.e. (E v n) > ksT s , and the rotational 
cooling, i.e. (E rot ) < kBT s . The final desorption an- 
gles are slightly larger than the experimental result [216|. 
The mean kinetic energy, however, is much larger for 
all ab initio molecular dynamics calculations than the 
experimental value for D2 of (Eki n ) exp = 0.17 eV > 
2k,BT s [190]. In the cluster calculations, molecules des- 
orbing from the isolated dihydride are even faster than 
molecules desorbing from the monohydride. It is remark- 
able that the results of the supercell and the cluster cal- 
culations are rather similiar for the mean final energies 
of hydrogen molecules desorbing from the monohydride 
state. 

The theoretical mean desorption angles are somewhat 
larger than the experimental values. However, due to 
the low number of trajectories obtained in the calcula- 
tions the discrepancies are not conclusive to rule out any 
particular mechanism. 

In conclusion, the analysis of the data of Tab. |h| gives 
no reason to disregard the prepairing mechanism with 
respect to any defect-mediated mechanism. According 
to the results of the ab initio molecular dynamics runs 
p9| , the isolated dihydride mechanism does even worse 
compared to the experiment than does the prepairing 
mechanism. And since the observed first-order kinetics 
in hydrogen desorption from Si(100) follows very natu- 
rally from the prepairing process, this still seems to be 
the most probable mechanism of desorption. Just re- 
cently there has been a very detailed temperature pro- 
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grammed desorption study of hydrogen desorption from 
Si(100) at low coverages which gave no indication that 
defect sites are involved in the desorption process [217| . 
Similar conclusions can be drawn from another recent 
study about hydrogen adsorption on Si(100) which al- 
lowed to discriminate between hydrogen adsorption on 
step and defect sites and on terrace sites |21S{ ]. 

The question still remains why the kinetic energy in 
desorption comes out much larger in the ab initio molec- 
ular dynamics runs than in the experiment. Uncertain- 
ties due to the bad statistics of the molecular dynamics 
runs or quantum dynamical effects are not believed to 
be responsible for the large difference f58,214|. In all the 
ab initio calculations the electrons are assumed to follow 
the motion of the nuclei adiabatically in their ground 
state. It has therefore been speculated that the excess 
energy at the transition state might be transferred to 
surface electronic excitations |5f||9l| thus taking energy 
away from the molecular degrees of freedom. Besides, the 
laser-induced desorption used in the experiment to deter- 
mine the translational energy in desorption fl9Cfl might 
not correspond to thermal desorption events described 
by the ab initio molecular dynamics runs. Recent DFT 
calculations indicate that laser melting of bulk silicon 
leads to a transition to a metallic liquid state 219 which 
could also occur during the laser-induced desorption used 
in the experiment. Furthermore, at surface tempera- 
tures above 900 K the clean Si(100) surface undergoes 
a semiconductor- metal transition p20j . At a metallic 
surface the barrier distribution might be different. In 
addition, the excitation of electron-hole pairs is much 
more probable leading to an effective dissipation channel 
for any excess energy at the transition state. Since all 
these considerations about the correct desorption mech- 
anism rely on one laser-induced desorption experiment, 
it is highly desirable that the mean kinetic energy of hy- 
drogen molecules desorbing from Si(100) is measured in- 
dependently by preferentially a thermal desorption ex- 
periment. 



B. Reactions of atomic hydrogen with 
hydrogen-passivated semiconductor surfaces 

The associative desorption reaction described in 
the previous section is of the so-called Langmuir- 
Hinshclwood (LH) type in which both rcactants are equi- 
librated with the substrate before the reaction. In the 
Eley-Rideal mechanism, on the other hand, an incoming 
gas phase species is assumed to react directly with an 
adsorbate and abstract it from the surface. In Fig. [l8| 
such an abstraction reaction is illustrated. Since in the 
ER reaction mechanism only one atom-surface bond has 
to be broken instead of two in the LH mechanism, the 
ER reactions show a significant exothermicity leading to 
large translational and internal energies in the desorbing 




FIG. 18. Example of a reactive trajectory for an H atom 
abstracting a D atom adsorbed on Si(100). The inci dent ki- 
netic energy of the H atom is 0.2 eV (from Ref. [ 229 1 ) . 



products. Although the concept has long been known, 
experimentally the Eley-Rideal mechanism has only re- 
cently been confirmed, at first for hydrogen abstraction 
reactions from metal surfaces [221-224]. In the mean- 
time also hydrogen abstraction reactions from Si surfaces 
have been observed [225,226]. These abstraction reac- 
tions have been found to be very efficient: the reaction 
probability for abstraction on the saturated monohydridc 
surface is 0.5±0.1 relative to the atomic adsorption prob- 
ability on the clean surface ]225[ |. Since the adsorption 
probability of hydrogen atoms on Si surfaces is of order 



unity [227], this means that roughly half of the incom- 
ing molecules pick up another hydrogen atom from the 
surface. However, only a small fraction of the incoming 
atoms directly hit the adsorbed hydrogen atoms. In or- 
der to explain the observed adsorption kinetics and the 
high abstraction probability, it was ther efore suggested 
that also indirect processes take place [225]. The im- 



pinging atom does not immediately pick up another atom 
from the surface but remains trapped at the surface while 
keeping most of its energy before the abstraction reac- 
tion. This mechanism, which has also been suggested for 
hydrogen abstraction re acti ons on metal surf aces [ 22S ] , is 
called "hot precursor" f225fl or "hot atom" |229[| mecha- 
nism and has first been addressed by Harris and Kasemo 



[230] 



Theoretically ER reactions have first been treated by 
low-dimensional quantum model studies in restricted ge- 
ometries [118,231,232]. These studies focused on the 
strong vibrational excitation of the molecules upon the 
abstraction reaction. Such a low-dimensional treatment 
does not allow the description of the hot atom process. 
In order to figure out the importance of these indirect 
processes, it is still unavoidable to rather perform clas- 
sical dynamical calculations with a high dimensionality 
instead of quantum calculations with a restricted dimen- 



26 



sionality. However, this is a reasonable approach since 
due to the large exothermicity of the reaction quantum 
effects should be negligible. Classical trajectory calcula- 
tions of the reaction dynamics of atomic hydrogen with 
the hydrogenated Si(100) surface have recently been per- 
formed by P. Kratzer ]229||62| ] . In these studies the PES 
has been described in a modified LEPS form. The en- 
ergy transfer to the Si substrate is taken into account by 
a collective mode in the surface oscillator model. In a 
first study, the LEPS parameters were derived from both 



experiment and ab initio calculations 1 229 ] . In a subse- 
quent treatment the PES parameters were all obtained 
by spin-polarized DFT-GGA calculations in the supercell 
approach 116211 . 



Figure |18| shows an example of a calculated reactive 
trajectory. An incoming H atom is directly abstracting 
an adsorbed D atom. The fast oscillations illustrate the 
strong vibrational excitation of the product HD molecule 
while the slow oscillations correspond to the rotational 
excitation. However, such a direct ER reaction is in fact 
less probable than an indirect hot atom abstraction re- 
action, as the ab initio molecular dynamics calculations 
reveal |62] ]. For H atoms incident with kinetic energies 
above 1 eV less than 20% of the produced HD molecules 
originate from direct ER reactions. 

Briefly summarizing the theoretical results, it is found 
that the abstraction probability, which agrees with the 
experiment within the experimental uncertainty, de- 
creases with increasing incident kinetic energy Most of 
the energy dissipated to the surface is transferred to vi- 
brations of the adsorbed hydrogens remaining on the sur- 
face since the energy transfer to the Si substrate atoms 
is rather low due to the large mass mismatch between 
H and Si. But the most important channel of energy 
disposal in the reaction is the kinetic energy of the hy- 
drogen molecules. Just recently also the H absorption in 
and adsorption on Cu(lll) has been studied by molecular 
dynamics calculations where the PES was derived from 
DFT calcula tion s and expanded in an effective-medium 
theory form | 233 |. These calculations confirm that ab ini- 
tio molecular dynamics studies are an indispensable tool 
to complement experimental information since contrary 
to the experiment the calculations allow a microscopic 
analysis of the reaction dynamics. 



V. LASER-INDUCED DESORPTION 

So far all reported ab initio dynamics studies employed 
the Born-Oppenheimer approximation, i.e., the electrons 
were assumed to follow the motion of the nuclei adia- 
batically. However, one very active field of experimen- 
tal research during the last years has been the study 
of desorption induced by electronic transitions, so-called 
DIET processes p34j. These processes can either be in- 
duced by electron impact (ESD) or by light (PSD). The 



abbreviations denote electron or photon stimulated des- 
orption, respectively. In particular, there have been nu- 
merous investigations of laser-induced desorption of small 
molecules from surfaces (for a review, see Ref. [235]). 

The theoretical description of the dynamics of desorp- 
tion processes involving electronically excited states is 
rather complex. In most studies the electron dynamics is 
not taken into account explicitly, instead the electronic 
transitions are treated in the Frank-Condon or sudden 



approximation [236,237]. In a quantum treatment of the 
dynamics of the nuclei this corresponds to the jump- 
ing of wave packets between different electronic states 
of the molecule-surface system. Usually only the elec- 
tronic ground state and one excited state are considered. 
In these studies the total energy is not conserved, i.e., it 
is assumed that the energy difference between the elec- 
tronic states is entirely taken away by the electrons. To 
overcome this problem, c ithe r the electron dynamics have 
to be treated explicitly [236 ] or the excit ation spectrum 
of the solid has to be taken into account [239]. 

Although the modelling of laser-induced desorption 
can be rather successful in elucidating dynamical aspects 



of the processes (see, e.g., Refs. [240 24I|), the potential 
curves for the electronic excited states had to be guessed 
since there were no accurate calculations for them. Fur- 
thermore, at metal surfaces the electronically excited 
states are strongly coupled to the continuous electronic 
excitation spectrum of the surface. It is questionable 
whether the description of a DIET process with just one 
sharp well-defined electronically excited state is reason- 
able in this strong-coupling regime. 

This situation is different for the excited states at in- 
sulator surfaces like oxides due to the large band gap. In 
addition, due to the strong ionicity of these systems a 
local description of the electronic structure is reasonable. 
This is important insofar as density functional theory is 
an electronic ground state theory and does not straight- 
forwardly allow the determination of excited state po- 
tentials (there are DFT calculations describing electronic 
excited states, though, see for example Ref. f92|). Hence 
one has to employ quantum chemical ab initio meth- 
ods which require the use of a cluster approach. In- 
deed such an approach has been followed recently for 
the first time. Kluner and coworkers have performed 
configuration interaction (CI) calculations for the con- 
struction of potential energy surfaces for electronically 
excited states describing the laser-induced desorption of 
NO from NiO(IOO) |||9]j24|] . In the CI calculations a 
NiOg~ cluster was embedded in a semi-infinite Madclung 
potential of ±2 point charges for the simulation of the 
NiO(f 00) surface. The calculated ground state potential 
did in fact not reproduce the experimentally found bind- 
ing energy of 0.52 eV |243|, confirming the well-known 
problems in determining binding energies at surfaces by 
cluster calculations pQ ). The ab initio minimum has 
therefore been scaled to fit the experimental data pi. 
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FIG. 19. Charge transfer PES of NO/NiO(100) as a func- 
tion of the NO center of mass distance R from the surface 
and the polar orientation of the molecule a. FC denotes the 
Franck-Condon point at which the wave packet propagation 
is started (from Ref. |B9|). 



justable parameter in the dynamics simulation. A value 
of r,. = 24 fs has been chosen which yields a desorption 
probability of 3.3% in agreement with typical experimen- 
tal data §jf|||]. 

Since the wave packet calculation are performed within 
a restricted dimensionality, quantitative agreement with 
the experiment [245] cannot be the ultimate goal of this 
theoretical study. More important is to gain a qualita- 
tive understanding of experimental trends. And indeed 
this study provides such a qualitative concept. The ex- 
cited state potential has a complicated shape as a func- 
tion of the polar orientation of the molecule. Two possi- 
ble pathways are sketched in Fig. [n] along which partial 
wave packets can propagate. This bifurcation of the wave 
packet causes a bimodality in the calculated velocity dis- 
tribution of desorbing molecules. This bimodality has 
also been found in the experiment [{245f . Hence the ab 
initio treatment of the laser-induced desorption relates 
qualitative trends in the experimental results to micro- 
scopic details of the potential energy surfaces. 



One specific charge transfer state in which one electron 
was transferred from the cluster to the NO molecule was 
used as a representative electronically excited state. The 
potential energy surfaces were determined as a function 
of the NO center of mass distance from the surface and 
the polar orientation of the molecule with respect to the 
surface normal. Figure |l9| shows the excited state poten- 
tial. 

The dynamics of the laser-induced desorption of NO 
from NiO(lOO) were simulated by the jumping wave 
packet method in which the molecular distance from 
the surface and the polar and azimuthal orientation 
of the molecule were considered explicitly. The laser- 
induced electronic excitation was modelled by putting the 
three-dimensional ground state wave function of the elec- 
tronic ground state onto the electronically excited PES. 
Thereby the center of the wave function is located at 
the Franck-Condon point FC shown on Fig. Il9[ The 
minimum of the excited state potential is closer to the 
surface than in the ground state potential leading to an 
acceleration of the wave packet towards the surface. This 
is called an Antoniewicz desorption scenario (244). The 
wave packet is then propagated for a certain residence 
time on the excited state potential before it is transferred 
in a second Franck-Condon transition back to the elec- 
tronic ground state. Molecules that have gained enough 
kinetic energy after the Franck-Condon transitions are 
able to desorb. The finite residence time is due to the 
interaction of the excited state with the electronic spec- 
trum of the substrate. The coupling strength of this 
interaction that causes the quenching of the electronic 
excitation is still unknown. Hence the mean residence 
time, the so-called resonance time r r , enters as an ad- 



VI. CONCLUSIONS AND OUTLOOK 

In this review ab initio dynamics studies of reactions 
on surfaces have been presented. In the last years the in- 
teraction between electronic structure calculations on the 
one side and dynamical calculations on the other side has 
been very fruitful. The availability of high-dimensional 
reliable potential energy surfaces has challenged the dy- 
namics community to improve their methods in order 
to perform high-dimensional dynamical studies on these 
potentials. Now quantum studies of the dissociation of 
hydrogen on surfaces are possible in which all six degrees 
of freedom of the molecule are treated dynamically. In 
this review I have tried to show that this achievement 
represents an important step forward in our understand- 
ing of the interaction of molecules with surfaces. Not 
only the quantitative agreement with experiment is im- 
proved, but also important qualitative concepts emerge 
from these high-dimensional calculations. 

However, what are the next steps to be taken in the 
application of ab initio dynamics studies? The six- 
dimensional dynamical treatment of hydrogen dissocia- 
tion represents the achievement of a goal that has long 
been pursued. But due to the unfavorable scaling of the 
quantum dynamical methods with the number of chan- 
nels there will probably be no seven- or eight-dimensional 
quantum calculations in the near future (although there 
are nevertheless some promising approaches [246]). Nei- 
ther will there soon be a six-dimensional quantum treat- 
ment of molecules heavier than hydrogen for the same 
reason. But for these heavier molecules quantum effects 
are not so important so that a classical treatment of the 
dynamics will be sufficient. 
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The important class of oxidation reactions, for exam- 
ple, will certainly soon be the subject of ab initio studies 
in which the dynamics is treated classically. In ab initio 
molecular dynamics simulation also the energy transfer 
to substrate phonons can be taken into account by ex- 
plicitly including the upper layers of the surface in the 
trajectory calculations. Still the quantum treatment of 
dissipation is an important goal because it would allow 
the energy transfer to phonons and to electron-hole pairs 
to be treated on an equal footing. This issue leads over to 
the treatment of electronically non-adiabatic processes. 
In this important field both the determination of excited 
state potentials as well as the description of non-adiabatic 
dynamics still represent a great challenge where much 
progress has to be made. 

Finally it should be mentioned that there is still room 
for low-dimensional model studies. In order to find out 
how the reaction dynamics depends on specific details of 
the PES it is usually very instructive to perform dynam- 
ical simulations in which certain features of the PES are 
changed in a controlled fashion within a restricted dimen- 
sionality. However, without any ab initio calculations it 
is often impossible to even guess the general shape of the 
PES for a more complicated reaction. This confirms the 
important role that ab initio dynamics calculations play 
for the advancement of our understanding of reactions on 
surfaces, even for model studies. 
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